Break apart the example below and save it under
the appropriate function names.  The first part
is a script and can be saved under any name you
wish, provided it has a .m suffix.


___________________


clear
echo on
%This program calculates the velocity and
%temperature distribution due to a heated 
%plane.  It uses the function fsolve to 
%obtain the two shooting parameters: the
%dimensionless shear stress and heat flux
%at the wall.  Because of the nature of the
%boundary layer equations, it is necessary to 
%have a very accurate initial guess for these
%parameters before the limit of integration
%etamax is made very large.  This is done by
%first choosing a small value of etamax and
%calculating the correct shooting parameters
%for this case, then using these as initial
%guesses for larger values of etamax.

format compact
global etamax pr umax

pr=1;
%It turns out that the initial value of etamax
%depends on the value of the Prandtl number
%that you use.  For large values of the Prandtl
%number we have to go to very small values of
%etamax initially to get stability, in part because
%the temperature profile decays over a very small
%region.  We thus adjust the values as follows:

n=max([round(pr^.5)-1,0]);
m=5;
etam=2.^[1-n:m];
%An initial guess that seems to work is:
finit=[1,-1]';

figure(1)
for i=1:n+m,
 etamax=etam(i)
 finit=fsolve('yinf',finit);
 shear_stress(i)=finit(1)
 heat_flux(i)=finit(2)
end;

%Finally, we look at the wall stress and flux
%as a function of etamax:
figure(2)
semilogx(etam,shear_stress,etam,heat_flux)
xlabel('eta limit')
ylabel('shear stress and heat flux')

%To conclude we replot the profiles over the
%range of eta = [1:6]:
etamax=6;
figure(1)
yinf(finit);

%And the maximum velocity is passed back through
%the variable umax:
umax

%This program took me about 25 minutes to write
%in its simple form (e.g., in order to get the 
%answer the first time) and another half an hour
%or so to comment and clean up.  Extending it to
%arbitrary large and small Prandtl numbers took
%another 45 minutes or so, most of it waiting for
%the program to execute.



___________________________


function y=yinf(finit)
global etamax umax

%We set our initial conditions for the integral:
deta=.05*min([1,etamax]);
eta=0;
f(1)=0;
f(2)=0;
f(3)=finit(1);
f(4)=1;
f(5)=finit(2);
i=1;

%We use a 4th order Runge-Kutta scheme from
%zero to some etamax.  This limit is passed
%in as a global parameter.
while eta < etamax,
 fplot(i,:)=[f(2),f(4)];
 etaplot(i)=eta;
 k1=deta*ydot(eta,f);
 k2=deta*ydot(eta+deta/2,f+k1/2);
 k3=deta*ydot(eta+deta/2,f+k2/2);
 k4=deta*ydot(eta+deta,f+k3);
 f=f+(k1+2*k2+2*k3+k4)/6;
 eta=eta+deta;
 i=i+1;
end;

%Some graphics to see how we are doing:
maxTV=max(fplot);
fplotscale=fplot*diag(1./maxTV);
plot(etaplot,fplotscale)
xlabel('eta')
ylabel('Scaled Velocity and Temperature')
drawnow

%We return the velocity and temperature at
%etamax to the function fsolve.  The two 
%shooting parameters will be adjusted until
%these output values are zero.
y(1)=f(2);
y(2)=f(4);

umax=max(fplot(:,1));



___________________________


function yp=ydot(eta,y)
%This function just gives the derivatives
%of the array of values.  The Prandtl
%number is passed in as a global variable.
global pr
yp(1)=y(2);
yp(2)=y(3);
yp(3)=-y(4)+2*y(2)^2-3*y(1)*y(3);
yp(4)=y(5);
yp(5)=-3*pr*y(1)*y(5);