Question

I need the matlab codes for following question

(1) (a). Solve the following second-order differential equations by a pair of first-order equations, xyʹʹ − yʹ − 8x³y³ = 0; with initial conditions y = 0.5 and yʹ = −0.5 at x = 1.

(b). Solve the problem in part (a) above using MATLAB built-in functions ode23 and ode45, within the range of 1 to 4, and compare with the exact solution of y = 1/(1 + x²)

[Hint: ode23 à 0.0456, ode45 à 0.0588]

(c). How can we improve the accuracy for the solutions obtained in parts (a) and (b) above?

Answer 1

%%Matlab code for solving ode
clear all
close all
%Answering question b.
%Initial conditions for ode
y0=[0.5;-0.5];
      
        %minimum and maximum x
        xspan=[1 4];
        %Solution of ODEs using ode45 matlab function
        sol1= ode45(@(x,y) odefcn1(x,y), xspan, y0);
        %Equally splitting x into small intervals
        x1 = linspace(xspan(1),xspan(end),501);
        %yy is the corresponding x y v and z
        yy1 = deval(sol1,x1);
        figure(1)
        hold on
        plot(x1,yy1(1,:),'linewidth',2)
      
        %minimum and maximum x
        xspan=[1 4];
        %Solution of ODEs using ode45 matlab function
        sol2= ode23(@(x,y) odefcn1(x,y), xspan, y0);
        %Equally splitting x into small intervals
        x1 = linspace(xspan(1),xspan(end),501);
        %yy is the corresponding x y v and z
        yy2 = deval(sol2,x1);
        figure(1)
     
        plot(x1,yy2(1,:),'linewidth',2)
     
      
        %exact solution
        y_ext=@(x) 1./(1+x.^2);
        plot(x1,y_ext(x1),'linewidth',2)
        hold off
      
title('Plot for y(x) vs. x solution')
xlabel('x')
ylabel('y(x)')
legend('ode45','ode23','Exact solution')
box on
grid on

fprintf('\tAt x=%2.f solution using ode45 is %2.5f.\n',x1(end),yy1(1,end))
fprintf('\n\tAt x=%2.f solution using ode23 is %2.5f.\n',x1(end),yy2(1,end))
fprintf('\n\tAt x=%2.f solution using exact analytic is %2.5f.\n',x1(end),y_ext(x1(end)))

%Function for evaluating the ODE
function dydx = odefcn1(x,y)

    eq1= y(2);
    eq2= (y(2)./x)+8.*x.^2.*(y(1)).^3;
  

    %Evaluate the ODE for our present problem
    dydx = [eq1;eq2];
  
end

%----------------------------End of Code-----------------------------%