Graph Equation

Solution

bc1.m

function res1=bc1(ya, yb)

% boundary conditions

res1=[ya(1); yb(1)-1];

end

bc2.m

function res2=bc2(ya, yb)

% boundary conditions

res2=[ya(1)-2; yb(1); yb(2)];

end

bc3.m

function res3=bc3(ya, yb)

% boundary conditions

res3=[ya(1); ya(2); yb(3); yb(4)-1];

end

ode1.m

function dydx1=ode1(x,y)

% y”+3yy’=0, y(0)=0, y(2)=1

% convert it into system of first order ODE

% let y1=y and y2=y’ then

% y1’=y2, y2’=-3*y1*y2

dydx1=[y(2); -3*y(1)*y(2)];

end

ode2.m

function dydx2=ode2(x,y)

% y”’=2y”+6xy, y(0)=2, y(5)=0, y'(5)=0

% convert it into system of first order ODE

% let y1=y, y2=y’, y3=y” then

% y1’=y2, y2’=y3, y3’=2*y3+6*x*y1

dydx2=[y(2); y(3); 2*y(3)+6*x*y(1)];

end

ode3.m

function dydx3=ode3(x,y)

% D4y=-4y^3/x, y(0)=0, y'(0)=0, y”(1)=0, y”'(1)=1

% convert it into system of first order ODE

% let y1=y, y2=y’, y3=y”, y4=y”’ then

% y1’=y2, y2’=y3, y3’=y4, y4’=-4y1^3/x

dydx3=[y(2); y(3); y(4); -4*y(1)^3/x];

end

solve.m

clc

%% Question 1

clear all

solinit = bvpinit(linspace(0,2,5),[0 0]);

sol = bvp4c(@ode1,@bc1,solinit);

x = linspace(0,2);

y = deval(sol,x);

figure

plot(x,y(1,:),’*’)

title(‘Solution 1’)

xlabel(‘x’);

ylabel(‘y’);

y = deval(sol,1.4);

fprintf(‘Solution 1\n’);

fprintf(‘y(x=1.4)=%f\n\n’,y(1));

fprintf(‘dy/dx(x=1.4)=%f\n\n’,y(2));

fprintf(‘d2y/dx2(x=1.4)=%f\n\n\n’,-3*y(1)*y(2));

%% Question 2

clear all

solinit = bvpinit(linspace(0,5,20),[0 0 0]);

sol = bvp4c(@ode2,@bc2,solinit);

x = linspace(0,5);

y = deval(sol,x);

figure

plot(x,y(1,:),’*’)

title(‘Solution 2’)

xlabel(‘x’);

ylabel(‘y’);

y = deval(sol,3);

fprintf(‘Solution 2\n’)

fprintf(‘y(x=3)=%f\n\n’,y(1));

fprintf(‘dy/dx(x=3)=%f\n\n’,y(2));

fprintf(‘d2y/dx2(x=3)=%f\n\n’,y(3));

fprintf(‘d3y/dx3(x=3)=%f\n\n\n’,2*y(3)+6*3*y(1));

%% Question 3

clear all

solinit = bvpinit(linspace(0.01,1,10),[1 1 1 1]);

sol = bvp4c(@ode3,@bc3,solinit);

x = linspace(0.01,1);

y = deval(sol,x);

figure

plot(x,y(1,:),’*’)

title(‘Solution 3’)

xlabel(‘x’);

ylabel(‘y’);

y = deval(sol,0.7);

fprintf(‘Solution 3\n’)

fprintf(‘y(x=0.7)=%f\n\n’,y(1));

fprintf(‘dy/dx(x=0.7)=%f\n\n’,y(2));

fprintf(‘d2y/dx2(x=0.7)=%f\n\n’,y(3));

fprintf(‘d3y/dx3(x=0.7)=%f\n\n’,y(4));