Setting up the initial conditions for BVP4C?

I'm trying to use BVP4C in a loop where the answer from the previous loop is the initial conditions for the next loop. Do I still use the Solinit function to get the initial conditions in the right format or should I do something else? Also, once the code gets to BVP4C, I get the following error:

Error using bvparguments (line 109) Error in calling BVP4C(ODEFUN,BCFUN,SOLINIT,P1,P2...): The derivative function ODEFUN should return a column vector of length 2.

I've seen some examples where the vector separator is a semicolon, but others where it appears to be a space. What exactly should I use?

The for loop that I used for testing is shown below:

%initial conditions
p(1) = pratio;
for i = 2:pnts
    p(i)= ((1/pnts/10)*i) + pratio;  %using x dimension
    Ma(i) = 0.0;
end

%start of full calculation loop, start with 2 iterations
%while (err>1*10^(-5))
for j = 1:2
ptot=0.0;
x=linspace(0,1,pnts);

%calculating q(x) eq. 26
for i = 1:pnts  
    Ao = 0.0;
    A_n =0.0;
    eta=acos(i/pnts);  
    cons = 1.0/sqrt((sinh(eps1)*sinh(eps1))+sin(eta)*sin(eta));
for k =1:200
      An = 0.0;
%to calculate An and Ao
%calculating the integral along fracture  
        for m = 1:(pnts-1)
        test =cos(2*acos(m/pnts));
        An = An + (1.0 - (p(m)*p(m)))*(cos(2*k*(acos(m/pnts)))/...
            sqrt(1.0 - m*m/pnts/pnts)); 
            if k==1
            Ao = Ao + (1.0 - (p(m)*p(m)))*(1.0 - m*m/pnts/pnts)^(-0.5);
            else
            end
        end  
  An = An*4.0/pi/(sinh(2.0*k*(eps_inf - eps1)));

      A_n = A_n + An;
      A_nn(k) = An;
  end
      Ao = (-2.0)/pi/(eps_inf - eps1)*Ao;

      q(i)= lamda*cons/p(i)*(A_n - Ao);

end

plot(x,q)

% For fixed well pressure:
%param(6)=pratio;            % if used for fixed well pressure condition

% How do I get q(x) in the right form?
solinit=bvpinit(linspace(0,1,pnts),[0.01; 0.9]); %0.01 intial Mach #, 0.9 initial pressure

%C3= (1.0 - gama*(y(1)^2.0));  % C3

sol=bvp4c(@frac,@fracbc,solinit);

%x=linspace(0,1,pnts);
y=deval(sol,x);

%store last values for error calculation
for i = 1:pnts
    p_1(i)= p(i);
    Ma_1(i) = Ma(i);
    ptot=ptot+p(i);
end

%write calculated values to arrays
for i = 1:pnts
    p(i)= y(2,i);
    Ma(i) = y(1,i);
end

%end of the while loop, for testing: a for loop
end

Here are the two nested functions that I forgot to put in:

   function res=fracbc(ya,yb) %#ok

    % y(1) is Mach #
    % y(2) is p, pressure

        Ao = 0.0;
        A_n =0.0;
    %calculating each the infinite series to 200 terms   
    for k =1:200
          An = 0.0;
    %to calculate An and Ao
    %calculating the integral along fracture  
            for m = 1:(pnts-1)
            An = An + (1.0 - (p(m)*p(m)))*(cos(2*k*(acos(m/pnts)))/...
                sqrt(1.0 - m*m/pnts/pnts)); 
                if k==1
                Ao = Ao + (1.0 - (p(m)*p(m)))*(1.0 - m*m/pnts/pnts)^(-0.5);
                else
                end
            end  
      An = An*4.0*k*(cosh(2.0*k*(eps_inf - eps1)))...
           /pi/(sinh(2.0*k*(eps_inf - eps1)));

        A_n = A_n + An;
        A_nn(k) = An;
    end
        Ao = (-2.0)/pi/(eps_inf - eps1)*Ao;
      

    M = (1.0/alpha/cfd/p(pnts))*(2.0*A_n - Ao);

    res = [ya(2)-pratio;yb(1)-M]; % This is for fixed pressure at well

    end

    %-------------------------------------------------------
    function dydx=frac(x,y)

    C3 = (1.0 - gama*(y(1).^2.0));
    C2 =sqrt(gama); 
    M_P = y(1)./y(2);

    B = (-y(1).*alpha - beta*y(2).*y(1).*y(1) - C2*q.*y(1).*y(2))/C3;

    A = (q./C2) - M_P.*B;

    dydx=[A;B];

    end

Answer 1

[Bill Greene]

From your call to bvpinit where you have this [0.01; 0.9] for the initial guess, it looks like you have 2 differential equations in your system. (In MATLAB, the ; means the entries that follow are in the next row so this is a column vector of length two.)

You didn't include the code for your frac function but apparently it isn't returning a column vector of length two, judging from the error message.

In MATLAB, there are two common ways to construct column vectors:

[1; 2; 3; 4] % create a column vector of length four directly

or

[1 2 3 4]' % create a row vector and then transpose it

Bill