How do I plot a graph from the following code

I need a graph of 'x vs y1' and 'x vs y0' as well as 'x vs y2', but the values of y0,y1,y2 has to be the iterated values.

% Fourth Runge Kutta for solving Blasius Equation

%-----------------------------------------------

%Blasius Equation : 2f"+ff'=0

%-----------------------------------------------

%Boundary Condition
%1. f'(0) = 0
%2. f'(inf) = 0
%3. f(0) = 0

%from shooting method
%4. f"(0) = 0.332

%-----------------------------------------------

%find: f,f" and f'''

%-----------------------------------------------

%Given

%eta = x

%f = y0,f'=y1,f"=y2,f'''= -(1/2)*y0*y2

%-----------------------------------------------

func1 = inline('y1','x','y0','y1','y2');

%-----------------------------------------------

func2 = inline('y2','x','y0','y1','y2');

%-----------------------------------------------

func3 = inline('-0.5*y0*y2','x','y0','y1','y2');

%-----------------------------------------------

%input: x = 0 , y0 = 0 , y1= 0

% y2 = 0.332 , total = 7 and h = 0.1

%-----------------------------------------------

x = input('\n Enter the value of x : ');

y0 = input( 'Enter the value of y0 : ');

y1 = input( 'Enter the value of y1 : ');

y2 = input( 'Enter the value of y2 : ');

total = input('Enter the value of total : ');

h = input( 'Enter the value of h : ');

fprintf('\n Solution with step size =%5.3f is:',h);

fprintf('\n x y0 y1 y2');

fprintf('\n%16.3f%16.3f%16.3f%16.3f',x,y0,y1,y2);

for i = 1:(total/h)

ak1y0 = func1(x,y0,y1,y2);

ak1y1 = func2(x,y0,y1,y2);

ak1y2 = func3(x,y0,y1,y2);

xx = x + h/2.;

yy0 = y0 + h*ak1y0/2.;

yy1 = y1 + h*ak1y1/2.;

yy2 = y2 + h*ak1y2/2.;

ak2y0 = func1(xx,yy0,yy1,yy2);

ak2y1 = func2(xx,yy0,yy1,yy2);

ak2y2 = func3(xx,yy0,yy1,yy2);

yy0 = y0 + h*ak2y0/2.;

yy1 = y1 + h*ak2y1/2.;

yy2 = y2 + h*ak2y2/2.;

ak3y0 = func1(xx,yy0,yy1,yy2);

ak3y1 = func2(xx,yy0,yy1,yy2);

ak3y2 = func3(xx,yy0,yy1,yy2);

xx = x + h;

yy0 = y0 + h*ak3y0;

yy1 = y1 + h*ak3y1;

yy2 = y2 + h*ak3y2;

ak4y0 = func1(xx,yy0,yy1,yy2);

ak4y1 = func2(xx,yy0,yy1,yy2);

ak4y2 = func3(xx,yy0,yy1,yy2);

y0 = y0 + (ak1y0 + 2.*ak2y0 + 2.*ak3y0 + ak4y0)*h/6.;

y1 = y1 + (ak1y1 + 2.*ak2y1 + 2.*ak3y1 + ak4y1)*h/6.;

y2 = y2 + (ak1y2 + 2.*ak2y2 + 2.*ak3y2 + ak4y2)*h/6.;

x = x + h;

fprintf('\n%16.3f%16.3f%16.3f%16.3f',x,y0,y1,y2);

end

Answer 1

[Walter Roberson]

all_x(i) = x;
all_y0(i) = y0;
all_y1(i) = y1;
all_y2(i) = y2;

Immediately before the x = x + h line.

Then after the loop,

plot(all_x, all_y0, 'k-', all_x, all_y1, 'b-', all_x, all_y2, 'g-')

[naygarp]

Thanks a lot. I got the graphs!

[naygarp]

Hey Walter, could you add another program to check the correct guess value i.e 'y2' by implementing linear interpolation or by secant method in this program itself. That would be very helpful.

1st guess, say y2= aplha1

2nd guess, say y2= alpha2

say phi(alpha1)= difference between the obtained boundary value and actual boundary value.

similarly,

say 'phi(alpha 2)'= difference between the obtained boundary value and the actual boundary value

after that

using Secant method for a better guess value 'alpha_n+1'

alpha_n+1 = alpha_n - (alpha_n - alpha_n-1)/(phi(alpha_n) - phi(alpha_n-1)*phi(alpha_n)

[Walter Roberson]

There have been a number of posts about secant method; no point my writing the code again.

[naygarp]

Hey Walter, how should I get the results rounded up to 5 decimal places in this particular code?

[Walter Roberson]

round() can be used to round to 5 decimal places, if round to nearest is acceptable instead of round-up. If you need round-up then

ceil(X*10^5)/10^5

As always, note that no binary floating point system can represent 0.1 exactly, just the same way that no fixed number of decimal places can represent 1/3 exactly. It is therefore not possible to round exactly to 5 decimal places, only to the nearest representable number.

[naygarp]

I meant to say that the values that appear in the command window after executing the code appears upto 3 decimal places like '1.235' . The modifications you made in the code

all_x(i) = x;

all_y0(i) = y0;

all_y1(i) = y1;

all_y2(i) = y2;

What should be the modification to get values upto 5 decimal places like '1.23567' if there is any?

[Walter Roberson]

There is no automatic method of displaying with 5 decimal places after the decimal point. You will need to use fprintf()

fmt = [repmat('%.5f ', 1, size(all_x,2)-1), '%.5f\n'];
fprintf(fmt, all_x.' );   %transpose is important

[naygarp]

Ok thanks I got it, if I use '%16.5f instead of '%16.3f' in printf(), then the results displayed would be upto 5 decimal places

[naygarp]

Hey Walter, could you help me in modifying the above code for another similar equation if you are familiar with using shooting method and runge-kutta 4th order numerical technique?