This is a boundary value problem (BVP), not an initial value problem (IVP). We can solve it using MATLAB's ODE45 solver by converting it from a BVP to an IVP through the shooting method. The general idea is the following:
- Disregard the boundary value(s) on the right. For this problem it means we ignore the condition that $f'(\infty) = 1.$
- Introduce new condition(s) on the left. In our case we've thrown out one condition on the right, so we have to introduce one condition on the left, say $f''(0) = 1$.
- Solve the initial value problem using the condition we introduced above. If we're very lucky the solution at the right hand endpoint agrees with the boundary conditions we discarded in step 1. More likely the solution will over- or undershoot the desired value.
- Make a sensible adjustment to the condition(s) we've introduced on the left and repeat from step 3. To decide what sort of adjustment is sensible you use your favourite root finding technique: bisection method, secant method, Newton's method etc.
Here are a few points (and some MATLAB code) for the implementation of the shooting method to this problem. Firstly, we need to rewrite our equation $ f''' + ff'' + \beta(1-f'^2) =0$ as a system of first order equations. This is not too difficult: $ f' = g, $ $ g' = h, $ $ h' = -fh - \beta(1-g^2).$ The corresponding function in MATLAB to be used with the ODE45 solver is something like:
function y_out = falkner_skan(x,y_in) global beta y_out(1,1) = y_in(2); y_out(2,1) = y_in(3); y_out(3,1) = -y_in(1)*y_in(3) - beta*(1-y_in(2)^2); end
And an implementation of the shooting method is:
clear global beta beta = 1; right_hand_endpoint = 5; stopping_distance = 0.01; shoot_parameter(1) = 1; initial_value = [0,0,shoot_parameter]; [x_out, y_out] = ode45(@falkner_skan,[0, right_hand_endpoint],initial_value); f_dash_endpoint(1) = y_out(end,2); shoot_parameter(2) = 1.1; n = 2; while abs(shoot_parameter(n) - shoot_parameter(n-1)) > stopping_distance initial_value = [0,0,shoot_parameter(n)]; [x_out, y_out] = ode45(@falkner_skan,[0, right_hand_endpoint],initial_value); f_dash_endpoint(n) = y_out(end,2); plot(x_out,y_out) pause shoot_parameter(n+1) = shoot_parameter(n) - (f_dash_endpoint(n)-1)*(shoot_parameter(n)-shoot_parameter(n-1))/(f_dash_endpoint(n)-f_dash_endpoint(n-1)); n = n+1; end
In the code above I've used the secant method to make my sensible guesses and instead of solving from $0$ to infinity I've chosen a value of $5$ as my right hand endpoint. The code is uncommented, but it's only 25 lines long and there is nothing mysterious going on in it so you should be able to follow it fairly easily. Let me know if you have questions.