%% Differential Equations
% This example shows how to use MATLAB® to formulate and solve several different types of differential equations. MATLAB offers several numerical algorithms to solve a wide variety of differential equations:
% Initial value problems
% Boundary value problems
% Delay differential equations
% Partial differential equations
%
%% Initial Value Problem
% vanderpoldemo is a function that defines the van der Pol equation
% \frac{d^2 y}{dt^2} - \mu \left( 1- y^2 \right) \frac{dy}{dt} + y = 0 .

type vanderpoldemo

%
% The equation is written as a system of two first-order ordinary differential equations (ODEs). These equations are evaluated for different values of the parameter \mu. For faster integration, you should choose an appropriate solver based on the value of \mu.
% For \mu =1, any of the MATLAB ODE solvers can solve the van der Pol equation efficiently. The ode45 solver is one such example. The equation is solved in the domain \left\lbrack 0,20\right\rbrack with the initial conditions y\left(0\right)=2 and \frac{\textrm{dy}}{\textrm{dt}}{\left|\right.}_{t=0} =0.

tspan = [0 20];
y0 = [2; 0];
Mu = 1;
ode = @(t,y) vanderpoldemo(t,y,Mu);
[t,y] = ode45(ode, tspan, y0);

% Plot solution
plot(t,y(:,1))
xlabel('t')
ylabel('solution y')
title('van der Pol Equation, \mu = 1')

%
% For larger magnitudes of \mu, the problem becomes stiff. This label is for problems that resist attempts to be evaluated with ordinary techniques. Instead, special numerical methods are needed for fast integration. The ode15s, ode23s, ode23t, and ode23tb functions can solve stiff problems efficiently.
% This solution to the van der Pol equation for \mu =1000 uses ode15s with the same initial conditions. You need to stretch out the time span drastically to \left\lbrack 0,3000\right\rbrack to be able to see the periodic movement of the solution.

tspan = [0, 3000];
y0 = [2; 0];
Mu = 1000;
ode = @(t,y) vanderpoldemo(t,y,Mu);
[t,y] = ode15s(ode, tspan, y0);

plot(t,y(:,1))
title('van der Pol Equation, \mu = 1000')
axis([0 3000 -3 3])
xlabel('t')
ylabel('solution y')

%
%% Boundary Value Problems
% bvp4c and bvp5c solve boundary value problems for ordinary differential equations.
% The example function twoode has a differential equation written as a system of two first-order ODEs. The differential equation is
% \frac{d^2 y}{dt^2} + |y|=0.

type twoode

%
% The function twobc has the boundary conditions for the problem: y\left(0\right)=0 and y\left(4\right)=-2.

type twobc

%
% Prior to calling bvp4c, you have to provide a guess for the solution you want represented at a mesh. The solver then adapts the mesh as it refines the solution.
% The bvpinit function assembles the initial guess in a form you can pass to the solver bvp4c. For a mesh of [0 1 2 3 4] and constant guesses of y\left(x\right)=1 and \textrm{y'}\left(x\right)=0, the call to bvpinit is:

solinit = bvpinit([0 1 2 3 4],[1; 0]);

%
% With this initial guess, you can solve the problem with bvp4c. Evaluate the solution returned by bvp4c at some points using deval, and then plot the resulting values.

sol = bvp4c(@twoode, @twobc, solinit);

xint = linspace(0, 4, 50);
yint = deval(sol, xint);
plot(xint, yint(1,:));
xlabel('x')
ylabel('solution y')
hold on

%
% This particular boundary value problem has exactly two solutions. You can obtain the other solution by changing the initial guesses to y\left(x\right)=-1 and \textrm{y'}\left(x\right)=0.

solinit = bvpinit([0 1 2 3 4],[-1; 0]);
sol = bvp4c(@twoode,@twobc,solinit);

xint = linspace(0,4,50);
yint = deval(sol,xint);
plot(xint,yint(1,:));
legend('Solution 1','Solution 2')
hold off

%
%% Delay Differential Equations
% dde23, ddesd, and ddensd solve delay differential equations with various delays. The examples ddex1, ddex2, ddex3, ddex4, and ddex5 form a mini tutorial on using these solvers.
% The ddex1 example shows how to solve the system of differential equations
% \begin{array}{l}
% y_1^{\prime } \left(t\right)=y_1 \left(t-1\right)\\
% y_2^{\prime } \left(t\right)=y_1 \left(t-1\right)+y_2 \left(t-0\ldotp 2\right)\\
% y_3^{\prime } \left(t\right)=y_2 \left(t\right)\ldotp
% \end{array}
% You can represent these equations with the anonymous function

ddex1fun = @(t,y,Z) [Z(1,1); Z(1,1)+Z(2,2); y(2)];

%
% The history of the problem (for t\le 0) is constant:
% \begin{array}{l}
% y_1 \left(t\right)=1\\
% y_2 \left(t\right)=1\\
% y_3 \left(t\right)=1\ldotp
% \end{array}
% You can represent the history as a vector of ones.

ddex1hist = ones(3,1);

%
% A two-element vector represents the delays in the system of equations.

lags = [1 0.2];

%
% Pass the function, delays, solution history, and interval of integration \left\lbrack 0,5\right\rbrack to the solver as inputs. The solver produces a continuous solution over the whole interval of integration that is suitable for plotting.

sol = dde23(ddex1fun, lags, ddex1hist, [0 5]);

plot(sol.x,sol.y);
title({'An example of Wille and Baker', 'DDE with Constant Delays'});
xlabel('time t');
ylabel('solution y');
legend('y_1','y_2','y_3','Location','NorthWest');

%
%% Partial Differential Equations
% pdepe solves partial differential equations in one space variable and time. The examples pdex1, pdex2, pdex3, pdex4, and pdex5 form a mini tutorial on using pdepe.
% This example problem uses the functions pdex1pde, pdex1ic, and pdex1bc.
%
% pdex1pde defines the differential equation
% \pi^2 \frac{\partial u}{\partial t} = \frac{\partial}{\partial x} \left(\frac{\partial u}{\partial x} \right).

type pdex1pde

%
% pdex1ic sets up the initial condition
% u(x,0) = \sin{\pi x}.

type pdex1ic

%
% pdex1bc sets up the boundary conditions
% u\left(0,t\right)=0,
% \pi e^{-t} +\frac{\partial }{\partial x}u\left(1,t\right)=0\ldotp

type pdex1bc

%
% pdepe requires the spatial discretization x and a vector of times t (at which you want a snapshot of the solution). Solve the problem using a mesh of 20 nodes and request the solution at five values of t. Extract and plot the first component of the solution.

x = linspace(0,1,20);
t = [0 0.5 1 1.5 2];
sol = pdepe(0,@pdex1pde,@pdex1ic,@pdex1bc,x,t);

u1 = sol(:,:,1);

surf(x,t,u1);
xlabel('x');
ylabel('t');
zlabel('u');

%
% Copyright 1984-2013 The MathWorks, Inc.