MATLAB: Loop possible? [Blasius Equation with fsolve]

blasiusboundary layerMATLAB

Hey ! So im pretty new to matlab, but I'm working my way into it. So now I have come across a problem, I'm not sure how to solve: I want to solve the blasius equation for a moving plate, thus different BC's than some might recognise (https://en.wikipedia.org/wiki/Blasius_boundary_layer) f''' + 1/2ff''=0, with BC's: f(0)=0, f'(0)=1 and f'(inf)=0. I know it's very easy to solve with the ODE (+shooting method) or rather the BVP command, but I need this example to apply it to a little more complex model, so I will need to do this with finite-differences(something similar to the keller box method). So since this will be boiling down to a problem formulation of a system of nonlinear equations, I figured why not use fsolve? I figured I would program a function containing all my BC's and equations, so I can just make a grid off "nodes" which will solve the equations at the corresponding node. However I seem to be stuck and or just stupid:
   %%Keller box Basius with fSolve
    eta0=0;
    etaInf=12;
    deltaEta=0.1; %stepsize
    N=(etaInf-eta0)/deltaEta; %number of nodes
    f=zeroes(N,1);g=ones(N,1);h=zeros(N,1); %freeing up space
    anfangswert=[0,1,0.3];
    %%call solver
    Sol = fsolve(blasiuskeller2(f,g,h,N,deltaEta), anfangswert);
    plot(eta0:edltaEta:etaInf, fval2);
  ---------------------------------------------------------------
    %%equations mit discretiszation 
    function fval = blasiuskeller2(f,g,h,N,deltaEta)
    % freeing up space
    fval1=zeros(N-1,1);fval2=zeros(N-1,1);fval3=zeros(N-1,1);
    % Define functions as F(X)=0
        for i=1:N-1
            fval1(i)=(f(i+1)-f(i))/deltaEta - (g(i+1)+g(i))/2;
            fval2(i)=(g(i+1)-g(i))/deltaEta - (h(i+1)+h(i))/2;
            fval3(i)=(h(i+1)-h(i))/deltaEta + (1/8)*(f(i+1)+f(i))*(h(i+1)+h(i));
        end;
    %BC's
    fval(:,1)=[f(1);g(1);fval2(1)]; %initial value
    for i=2:(N-1)
        fval(:,i)=[fval1(i-1);fval3(i-1);fval2(i)];
    end;
    fval(:,N)=[fval1(N-1);1-g(N);fval2(N-1)]; %Boundary value
Pretty sure this is possible, but not sure if I'm on the right track here. Any help would be greatly appreciated!

Best Answer

I did not check your equations, but already your call to fsolve is incorrect:
anfangswert=horzcat(zeros(1,1:N),ones(1,1:N),0.3*ones(1,1:N));
sol = fsolve(@(x)blasiuskeller2(x(1:N),x(N+1:2*N),x(2*N+1:3*N),N,deltaEta), anfangswert);
Best wishes
Torsten.
Related Question