function [ret] = hydraulic_system(request,z,za,zb,zc,t,p,lambda,lambda_p)
% Problem definition for parameter dependent problem in Manual
global G1 kr kx ha pwp PB LB rc r0 T Mb L Psis

switch request
    case 'n'
        ret = 3;
    case 'orders'
        ret = [ 2 0 1];
    case 'problem'
         ret =   [
                  (z(1,3) - 2*pi*r0*kr/kx*(z(1,1) -  ha/((ha/pwp)^PB + z(2,1)/LB)^(1/PB))) ; 
                  (z(2,1) - (p(1) - G1*z(1,3)*kx));   % M0  p(1)= Mrc
                  (z(3,2) - (z(2,1) +  z(1,3)*kx*(2*G1-1)/(pi*(rc^2-r0^2)))) % mean MFP 
                  ];   
                
% syms  z1 z2 ha pwp x21 LB PB r0 kr kx pi
% diff(z1 - 2*pi*r0*kr/kx*(z2 -  ha/((ha/pwp)^PB + x21/LB)^(1/PB)),x21,1)
             
    case 'jacobian'
        %DON'T CHANGE THIS LINE:
        ret = zeros(length(feval(mfilename,'problem',zeros(feval(mfilename,'n'),max(feval(mfilename,'orders'))+1),[],[],[],0,zeros(feval(mfilename,'parameters'),1))),feval(mfilename,'n'),max(feval(mfilename,'orders'))+1);
        ret(1, 1, 1) = -2*pi*r0*kr/kx;
        ret(1, 1, 3) = 1;
        ret(1, 2, 1) = -(2*ha*kr*pi*r0)/(LB*PB*kx*((ha/pwp)^PB + z(2,1)/LB)^(1/PB + 1));
        ret(2, 1, 3) = G1*kx;
        ret(2, 2, 1) = 1;
        ret(3, 1, 3) = kx*(2*G1-1)/(pi*(rc^2-r0^2));
        ret(3, 2, 1) = -1;
        ret(3, 3, 2) = 1;
        if sum(isnan(ret))+sum(isinf(ret))>0
            disp('asd')
        end
    case 'interval'
        ret = [0, L];
    case 'linear'
        ret = 0;
    case 'parameters'
        ret = 2;
    case 'c'
        ret = [];
    case 'BV'
        ret=[
            za(1,2);
            zb(1,1) - p(2); % p(2) = PX;
            za(3,1);  
            zb(1,2)*kx - T; 
            zb(3,1)/L -  Mb         
            ];

    case 'dBV'
        %DON'T CHANGE THIS LINE:
        ret = zeros(max(length(feval(mfilename,'c')),2-length(feval(mfilename,'c'))),length(feval(mfilename,'problem',zeros(feval(mfilename,'n'),max(feval(mfilename,'orders'))+1),[],[],[],0,zeros(feval(mfilename,'parameters'),1))),feval(mfilename,'n'),max(feval(mfilename,'orders')));
        ret(1,1,1,2) = 1;
        ret(2,2,1,1) = 1;
        ret(1,3,3,1) = 1;
        ret(2,4,1,2) = kx;
        ret(2,5,3,1) = 1/L;
    case 'dP'
        ret = [              
                  0,  0  ;
                  -1, 0  ;
                  0,  0
               ];  
    case 'dP_BV'
        ret =   [ 
        0 , 0; 
        0 , -1;
        0 ,0 ; 
        0, 0 ; 
        0, 0 
        ] ;       
    case 'initProfile'
        ret.initialMesh    = linspace(0,L,L*6);
        ret.initialValues  = [ones(1,L*6)*Psis; ones(1,L*6)*Mb; zeros(1,L*6)];             
        ret.parameters     = [Mb*1.1, Psis*1.5] ;
    case 'EVP'
        ret = 0;
    case 'dLambda'
        ret = 0;
    case 'pathfollowing'
        ret.activate=0;
    otherwise
        ret = 0;                
end

end