Solve nonstiff differential equations — medium order method (2024)

Simple ODEs that have a single solution component can be specified as an anonymous function in the call to the solver. The anonymous function must accept two inputs (t,y), even if one of the inputs is not used in the function.

Solve the ODE

$$y^{\prime }=2t.$$

Specify a time interval of [0 5] and the initial condition y0 = 0.

tspan = [0 5];y0 = 0;[t,y] = ode45(@(t,y) 2*t, tspan, y0);

Plot the solution.

plot(t,y,'-o')

The van der Pol equation is a second-order ODE

$$y_1^{\prime \prime }-\mu (1-y_1^2)y_1^{\prime }+y_1=0,$$

where $$\mu >0$$ is a scalar parameter. Rewrite this equation as a system of first-order ODEs by making the substitution $$y_1^{\prime }=y_2$$. The resulting system of first-order ODEs is

$$\begin{array}{cl}{y}_1^{\prime }& =y_2\\ y_2^{\prime }& =\mu (1-y_1^2)y_2-y_1.\end{array}$$

The function file vdp1.m represents the van der Pol equation using $$\mu =1$$. The variables $$y_1$$ and $$y_2$$ are the entries y(1) and y(2) of a two-element vector dydt.

function dydt = vdp1(t,y)%VDP1 Evaluate the van der Pol ODEs for mu = 1%% See also ODE113, ODE23, ODE45.% Jacek Kierzenka and Lawrence F. Shampine% Copyright 1984-2014 The MathWorks, Inc.dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];

Solve the ODE using the ode45 function on the time interval [0 20] with initial values [2 0]. The resulting output is a column vector of time points t and a solution array y. Each row in y corresponds to a time returned in the corresponding row of t. The first column of y corresponds to $$y_1$$, and the second column corresponds to $$y_2$$.

[t,y] = ode45(@vdp1,[0 20],[2; 0]);

Plot the solutions for $$y_1$$ and $$y_2$$ against t.

plot(t,y(:,1),'-o',t,y(:,2),'-o')title('Solution of van der Pol Equation (\mu = 1) with ODE45');xlabel('Time t');ylabel('Solution y');legend('y_1','y_2')

ode45 works only with functions that use two input arguments, t and y. However, you can pass extra parameters by defining them outside the function and passing them in when you specify the function handle.

Solve the ODE

$$y^{\prime \prime }=\frac{A}{B}ty.$$

Rewriting the equation as a first-order system yields

$$\begin{array}{cl}{y}_1^{\prime }& =y_2\\ y_2^{\prime }& =\frac{A}{B}t{y}_1.\end{array}$$

odefcn, a local function included at the end of this example, represents this system of equations as a function that accepts four input arguments: t, y, A, and B.

function dydt = odefcn(t,y,A,B) dydt = zeros(2,1); dydt(1) = y(2); dydt(2) = (A/B)*t.*y(1);end

Solve the ODE using ode45. Specify the function handle so that it passes the predefined values for A and B to odefcn.

A = 1;B = 2;tspan = [0 5];y0 = [0 0.01];[t,y] = ode45(@(t,y) odefcn(t,y,A,B), tspan, y0);

Plot the results.

plot(t,y(:,1),'-o',t,y(:,2),'-.')

function dydt = odefcn(t,y,A,B) dydt = zeros(2,1); dydt(1) = y(2); dydt(2) = (A/B)*t.*y(1);end

For simple ODE systems with one equation, you can specify y0 as a vector containing multiple initial conditions. This technique creates a system of independent equations through scalar expansion, one for each initial value, and ode45 solves the system to produce results for each initial value.

Create an anonymous function to represent the equation $\mathit{f}\left(\mathit{t},\mathit{y}\right)=-2\mathit{y}+2\cos \left(\mathit{t}\right)\sin \left(2\mathit{t}\right)$. The function must accept two inputs for t and y.

yprime = @(t,y) -2*y + 2*cos(t).*sin(2*t);

Create a vector of different initial conditions in the range $\left[-5,5\right]$.

y0 = -5:5; 

Solve the equation for each initial condition over the time interval $\left[0,3\right]$ using ode45.

tspan = [0 3];[t,y] = ode45(yprime,tspan,y0);

Plot the results.

plot(t,y)grid onxlabel('t')ylabel('y')title('Solutions of y'' = -2y + 2 cos(t) sin(2t), y(0) = -5,-4,...,4,5','interpreter','latex')

This technique is useful for solving simple ODEs with several initial conditions. However, the technique also has some tradeoffs:

For more information on this technique, see Solve System of ODEs with Multiple Initial Conditions.

The van der Pol equation is a second order ODE

$$y_1^{\prime \prime }-\mu (1-y_1^2)y_1^{\prime }+y_1=0.$$

Solve the van der Pol equation with $$\mu =1$$ using ode45. The function vdp1.m ships with MATLAB® and encodes the equations. Specify a single output to return a structure containing information about the solution, such as the solver and evaluation points.

tspan = [0 20];y0 = [2 0];sol = ode45(@vdp1,tspan,y0)

sol = struct with fields: solver: 'ode45' extdata: [1x1 struct] x: [0 1.0048e-04 6.0285e-04 0.0031 0.0157 0.0785 0.2844 0.5407 0.8788 1.4032 1.8905 2.3778 2.7795 3.1285 3.4093 3.6657 3.9275 4.2944 4.9013 5.3506 5.7998 6.2075 6.5387 6.7519 6.9652 7.2247 7.5719 8.1226 8.6122 9.1017 9.5054 ... ] (1x60 double) y: [2x60 double] stats: [1x1 struct] idata: [1x1 struct]

Use linspace to generate 250 points in the interval [0 20]. Evaluate the solution at these points using deval.

x = linspace(0,20,250);y = deval(sol,x);

Plot the first component of the solution.

plot(x,y(1,:))

Extend the solution to $$t_f=35$$ using odextend and add the result to the original plot.

sol_new = odextend(sol,@vdp1,35);x = linspace(20,35,350);y = deval(sol_new,x);hold onplot(x,y(1,:),'r')