6  MATLAB’s In-Built Procedures

So far, the three main iterative methods have been developed that solve IVPs numerically. MATLAB, however, has its own built-in procedures that can solve IVPs with a combination of several methods. The two main ones are ode23 (which uses a combination of a second and third order RK methods) and ode45 (which uses a combination of a fourth and fifth order RK methods).

Both ode45 and ode23 are hybrid methods and use adaptive meshing, this means that the time span grid is not necessarily uniform, but it changes depending on the gradients; if the gradient is large at some point, then the stepsize will be small to capture these drastic changes.

The following MATLAB code solves the following set of IVPs on the interval \([0,1]\) using ode45: \[\begin{align*} & \frac{\mathrm{d} u}{\mathrm{d} t}=2u+v+w+\cos(t), & \quad u(0)=0 \\ & \frac{\mathrm{d} v}{\mathrm{d} t}=\sin(u)+\mathrm{e}^{-v+w}, & \quad v(0)=1 \\ & \frac{\mathrm{d} w}{\mathrm{d} t}=uv-w, & \quad w(0)=0. \end{align*}\]

function IVP_InBuilt

%% Solve a set of first order IVPs using In-Built codes

% This code solves a set of IVP when written explicitly
% on the interval [t0,tf] subject to the initial conditions
% y(0)=y0.  The output will be the graph of the solution(s)
% and the vector value at the final point tf.  Note that the
% IVPs do not need to be linear or homogeneous.

%% Lines to change:

% Line 28   : t0 - Start time
% Line 31   : tf - End time
% Line 43   : T_Span - Time span for evaluation
% Line 46   : y0 - Vector of initial values
% Line 86+  : Which functions to plot, remembering to assign
%             a colour, texture and legend label
% Line 100+ : Set of differential equations written
%             explicitly.  These can also be non-linear and
%             include forcing terms.  These equations can
%             also be written in matrix form if the
%             equations are linear.

%% Set up input values

% Start time
t0=0;

% End time
tf=1;

% Time span
% In-built methods tend to use adaptive meshing; decreasing
% the stepsize near locations with drastic derivative
% changes and increasing near small derivative changes.
% Sometimes this is not desired but a uniform meshing is
% requiredfrom the start time t0 to the end time tf being
% split into N equal sub intervals.  This can be changed
% here:
% Adaptive meshing:  T_Span=[t0 tf]
% Specific meshing:  T_Span=linspace(t0,tf,N)
T_Span=[t0 tf];

% Column vector initial values y0=y(t0)
y0=[0;1;0];

%% Set up IVP solver parameters

% Number of differential equations
K=length(y0);

%% Use solver

% Set the solver tolerance
tol=odeset('RelTol',1e-6);

% Solve the IVP using ode45 or ode23
[T,Y]=ode45(@(t,y) DYDT(t,y,K),T_Span,y0,tol);

% Convert T and Y to columns for consistency
T=T';
Y=Y';

%% Setting plot parameters

% Clear figure
clf

% Hold so more than one line can be drawn
hold on

% Turn on grid
grid on

% Setting font size and style
set(gca,'FontSize',20,'FontName','Times')
set(legend,'Interpreter','Latex')

% Label the axes
xlabel('$t$','Interpreter','Latex')
ylabel('$\mathbf{y}(t)$','Interpreter','Latex')

% Plot the desried solutions.  If all the solutions are
% needed, then consider using a for loop in that case
plot(T,Y(1,:),'-b','LineWidth',2,'DisplayName','$y_1(t)$')
plot(T,Y(2,:),'-r','LineWidth',2,'DisplayName','$y_2(t)$')
plot(T,Y(3,:),'-k','LineWidth',2,'DisplayName','$y_3(t)$')

% Display the values of the vector y at tf
disp(strcat('The vector y at t=',num2str(tf),' is:'))
disp(Y(:,end))

end

function [dydt]=DYDT(t,y,K)

% When the equation are written in explicit form

dydt=zeros(K,1);

dydt(1)=2*y(1)+y(2)+y(3)+cos(t);

dydt(2)=sin(y(1))+exp(-y(2)+y(3));

dydt(3)=y(1)*y(2)-y(3);

% If the set of equations is linear, then these can be
% written in matrix form as dydt=A*y+b(t).  For example, if
% the set of equations is:
% dudt = 7u - 2v +  w + exp(t)
% dvdt = 2u + 3v - 9w + cos(t)
% dwdt =      2v + 5w + 2
% Then:
% A=[7,-2,1;2,3,-9;0,2,5];
% b=@(t) [exp(t);cos(t);2];
% dydt=A*y+b(t)

end