{"id":1589,"date":"2023-10-03T06:23:36","date_gmt":"2023-10-03T10:23:36","guid":{"rendered":"https:\/\/blogs.mathworks.com\/matlab\/?p=1589"},"modified":"2023-10-05T19:15:02","modified_gmt":"2023-10-05T23:15:02","slug":"the-new-solution-framework-for-ordinary-differential-equations-odes-in-matlab-r2023b","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/matlab\/2023\/10\/03\/the-new-solution-framework-for-ordinary-differential-equations-odes-in-matlab-r2023b\/","title":{"rendered":"The new solution framework for Ordinary Differential Equations (ODEs) in MATLAB R2023b"},"content":{"rendered":"
Along with <\/span>linear algebra<\/span><\/a>, one of the iconic features of MATLAB in my mind is how it handles <\/span>ordinary differential equations (ODEs)<\/span><\/a>. ODEs have been part of MATLAB almost since the very beginning.<\/span><\/div>
One of the features of how MATLAB traditionally allows users to solve ODEs is that it provides a suite of functions. For many years, there were 7 routines in the MATLAB suite, described in <\/span>this 2014 post by Cleve Moler<\/span><\/a> and in more detail in the 1997 paper, <\/span>The MATLAB ODE Suite<\/span><\/a>. In 2021b, 2 new high-order methods were added to the suite, <\/span>ode78<\/span><\/a> and <\/span>ode89<\/span><\/a> based on the algorithms described in the paper <\/span>Numerically optimal Runge\u2013Kutta pairs with interpolants<\/span><\/a>. <\/span><\/div>
The design of the suite was elegant and has served the community well for over 25 years! However, MATLAB has evolved a lot since it was designed and user-expectations have evolved with it. It was felt that it was time for a fresh look at how to solve ODEs in MATLAB, one that would additionally support our future plans in a modern and elegant manner. <\/span><\/div>
Before I dive into the new interface, I'd like to make it clear that <\/span>the existing suite of ODE functions are not going anywhere! <\/span>Millions of lines of code make use of <\/span>ode45<\/span><\/a> and friends and we have no plans on doing anything that would break that code. This is about adding functionality, not taking anything away. <\/span><\/div>
Furthermore, the focus of this release is the new interface. There are no new solvers or any fundamentally new functionality....yet! However, some things will be significantly easier to do using the new interface.<\/span><\/div>
With that said, let's take a look at how the new design looks.<\/span><\/div>

Solving ODEs in MATLAB the OOP way<\/span><\/h2>
Say we want to solve the ODE<\/span><\/div>
y<\/mi><\/mrow>\u2032<\/mo><\/mrow><\/msup>=<\/mo>2<\/mn>t<\/mi><\/mrow><\/math><\/span><\/div>
with the initial condition <\/span>y<\/mi><\/mrow>0<\/mn><\/mrow><\/msub>=<\/mo>0<\/mn><\/mrow><\/math><\/span> over the interval <\/span>t<\/mi>=<\/mo>[<\/mo>0<\/mn>,<\/mo>10<\/mn><\/mrow>]<\/mo><\/mrow><\/mrow><\/math><\/span><\/div>
With the new interface, we can set up and solve this problem as follows<\/span><\/div>
F = ode(ODEFcn=@(t,y) 2*t,InitialTime=0,InitialValue=0); <\/span>% Set up the problem by creating an ode object<\/span><\/span><\/div><\/div>
sol = solve(F,0,10); <\/span>% Solve it over the interval [0,10]<\/span><\/span><\/div><\/div>
plot(sol.Time,sol.Solution,<\/span>\"-o\"<\/span>) <\/span>% Plot the solution<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
An alternative way to proceed would have been to start with an empty <\/span>ode object<\/span><\/a> and add one property at a time:<\/span><\/div>
F = ode; <\/span>% Empty ode object called F<\/span><\/span><\/div><\/div>
F.ODEFcn = @(t,y) 2*t; <\/span>% Add the function we want to solve to F<\/span><\/span><\/div><\/div>
F.InitialTime = 0; <\/span>% Add the initial time to F<\/span><\/span><\/div><\/div>
F.InitialValue = 0; <\/span>% Add the initial value to F<\/span><\/span><\/div><\/div><\/div>
Once our problem is set up, we pass it to the <\/span>solve<\/span> function.<\/span><\/div>
sol = solve(F,0,10); <\/span><\/span><\/div><\/div><\/div>
The solution is an <\/span>ODEResults<\/span><\/a> object with two properties: <\/span>Time<\/span> and <\/span>Solution<\/span><\/div>
sol<\/span><\/span><\/div>
sol = <\/span><\/div>
ODEResults with properties:\n\n Time: [0 0.2500 0.5000 0.7500 1 1.2500 1.5000 1.7500 2 2.2500 2.5000 2.7500 3 3.2500 3.5000 3.7500 4 4.2500 4.5000 4.7500 5 5.2500 5.5000 5.7500 6 6.2500 6.5000 6.7500 7 7.2500 7.5000 7.7500 8 8.2500 8.5000 8.7500 9 9.2500 9.5000 9.7500 10]\n Solution: [0 0.0625 0.2500 0.5625 1.0000 1.5625 2.2500 3.0625 4 5.0625 6.2500 7.5625 9 10.5625 12.2500 14.0625 16 18.0625 20.2500 22.5625 25 27.5625 30.2500 33.0625 36 39.0625 42.2500 45.5625 49 52.5625 56.2500 60.0625 64 68.0625 \u2026 ] (1\u00d741 double)\n<\/div><\/div><\/div><\/div><\/div><\/div>
You access either property with . notation<\/span><\/div>
sol.Solution<\/span><\/span><\/div>
ans = <\/span>1\u00d741<\/span><\/div><\/div>
0 0.0625 0.2500 0.5625 1.0000 1.5625 2.2500 3.0625 4.0000 5.0625 6.2500 7.5625 9.0000 10.5625 12.2500 14.0625 16.0000 18.0625 20.2500 22.5625 25.0000 27.5625 30.2500 33.0625 36.0000 39.0625 42.2500 45.5625 49.0000 52.5625 56.2500 60.0625 64.0000 68.0625 72.2500 76.5625 81.0000 85.5625 90.2500 95.0625 100.0000\n<\/div>
<\/div>
<\/div><\/div><\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div>
sol.Time<\/span><\/span><\/div>
ans = <\/span>1\u00d741<\/span><\/div><\/div>
0 0.2500 0.5000 0.7500 1.0000 1.2500 1.5000 1.7500 2.0000 2.2500 2.5000 2.7500 3.0000 3.2500 3.5000 3.7500 4.0000 4.2500 4.5000 4.7500 5.0000 5.2500 5.5000 5.7500 6.0000 6.2500 6.5000 6.7500 7.0000 7.2500 7.5000 7.7500 8.0000 8.2500 8.5000 8.7500 9.0000 9.2500 9.5000 9.7500 10.0000\n<\/div>
<\/div>
<\/div><\/div><\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>

Automatic solver selection<\/span><\/h2>
One of the first things to note in the above workflow is that we didn't choose a solver. That is, I didn't have to think about if I should use <\/span>ode45<\/span>, <\/span>ode78<\/span> etc. All I did was state the problem mathematically and ask MATLAB to solve it. Let's take a closer look at the details.<\/span><\/div>
I set up my problem like this:<\/span><\/div>
F = ode(ODEFcn=@(t,y) 2*t,InitialTime=0,InitialValue=0);<\/span><\/span><\/div><\/div><\/div>
F is an <\/span>ode<\/span> object:<\/span><\/div>
class(F)<\/span><\/span><\/div>
ans = <\/span>'ode'<\/div><\/div><\/div><\/div><\/div><\/div>
I can see the details by evaluating F<\/span><\/div>
F<\/span><\/span><\/div>
F = <\/span><\/div>
ode with properties:\n\n Problem definition<\/strong>\n ODEFcn: @(t,y)2*t\n InitialTime: 0\n InitialValue: 0\n\n Solver properties<\/strong>\n AbsoluteTolerance: 1.0000e-06\n RelativeTolerance: 1.0000e-03\n Solver: auto\n SelectedSolver: ode45\n\n Show all properties\n\n<\/div><\/div><\/div><\/div><\/div><\/div>
The object display separates the mathematical problem definition from how we are going to solve it. By default, the <\/span>Solver<\/span> property of the class is set to <\/span>\"auto\"<\/span> which means that MATLAB attempts to choose a suitable solver based on various properties of the problem we've asked it to solve. In this case it has chosen <\/span>ode45<\/span> which is a pretty safe bet for many problems. <\/span><\/div>
Ask for a tighter tolerance for this problem, however, and it switches to using <\/span>ode78<\/span>.<\/span><\/div>
F.RelativeTolerance = 1e-7<\/span><\/span><\/div>
F = <\/span><\/div>
ode with properties:\n\n Problem definition<\/strong>\n ODEFcn: @(t,y)2*t\n InitialTime: 0\n InitialValue: 0\n\n Solver properties<\/strong>\n AbsoluteTolerance: 1.0000e-06\n RelativeTolerance: 1.0000e-07\n Solver: auto\n SelectedSolver: ode78\n\n Show all properties\n\n<\/div><\/div><\/div><\/div><\/div><\/div>
At this stage, MATLAB hasn't solved anything. We've just set up our problem and made decisions about <\/span>how<\/span> we are going to solve it. We need the <\/span>solve<\/span> function to complete the job. <\/span><\/div>
One thing to bear in mind is that the solver that <\/span>\"auto\"<\/span> chooses may change in future releases for a number of reasons. For example, we may improve the heuristics used to select the best solver or maybe add new solvers that do a better job for a given problem. <\/span><\/div>
If you want to override what MATLAB chooses and fix the solver type you can do that as follows<\/span><\/div>
F.Solver = <\/span>\"ode45\"<\/span> <\/span>% Force the framework to use ode45. <\/span><\/span><\/div>
F = <\/span><\/div>
ode with properties:\n\n Problem definition<\/strong>\n ODEFcn: @(t,y)2*t\n InitialTime: 0\n InitialValue: 0\n\n Solver properties<\/strong>\n AbsoluteTolerance: 1.0000e-06\n RelativeTolerance: 1.0000e-07\n Solver: ode45\n\n Show all properties\n\n<\/div><\/div><\/div><\/div><\/div><\/div>

Trying a stiff ODE problem<\/span><\/h2>
I was curious about how it would handle a stiff ODE and so fed it the example given by Cleve in his blog post <\/span>Ordinary Differential Equations, Stiffness<\/span><\/a>. <\/span><\/div>
delta = 0.0001;<\/span><\/span><\/div><\/div>
CleveF = ode(ODEFcn=@(t,y) y^2 - y^3);<\/span><\/span><\/div><\/div>
CleveF.InitialTime = 0;<\/span><\/span><\/div><\/div>
CleveF.InitialValue = delta;<\/span><\/span><\/div><\/div>
sol = solve(CleveF,0,2\/delta);<\/span><\/span><\/div><\/div>
plot(sol.Time,sol.Solution,<\/span>'-o'<\/span>)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
Looks reasonable. Following Cleve's suggestion, we zoom in on the steady state that begins at <\/span>x=1 <\/span>and see that the solver is working hard to do its job, just as it did when Cleve explicitly chose <\/span>ode45<\/span>. <\/span><\/div>
plot(sol.Time,sol.Solution,<\/span>'-o'<\/span>)<\/span><\/span><\/div><\/div>
xlim([10006.3 10112.1])<\/span><\/span><\/div><\/div>
ylim([0.99615 1.00250])<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
Sure enough, the <\/span>\"auto\" <\/span>option has also chosen ode45 in this case. <\/span><\/div>
CleveF.SelectedSolver<\/span><\/span><\/div>
ans = <\/span><\/div>
SolverID enumeration\n\n ode45\n<\/div><\/div><\/div><\/div><\/div><\/div>
Using the old solver suite, I'd now have to dig into the documentation and read about the stiff solvers available to me before trying those out. With the new interface, however, I just tell it that I think the problem is stiff and it will choose a stiff solver for me. <\/span><\/div>
CleveF.Solver = <\/span>\"stiff\"<\/span><\/span><\/div>
CleveF = <\/span><\/div>
ode with properties:\n\n Problem definition<\/strong>\n ODEFcn: @(t,y)y^2-y^3\n InitialTime: 0\n InitialValue: 1.0000e-04\n Jacobian: []\n\n Solver properties<\/strong>\n AbsoluteTolerance: 1.0000e-06\n RelativeTolerance: 1.0000e-03\n Solver: stiff\n SelectedSolver: ode15s\n\n Show all properties\n\n<\/div><\/div><\/div><\/div><\/div><\/div>
This time it has selected <\/span>ode15s<\/span><\/a> for me. Let's see how that does <\/span><\/div>
stiffSol = solve(CleveF,0,2\/delta);<\/span><\/span><\/div><\/div>
plot(stiffSol.Time,stiffSol.Solution,<\/span>'-o'<\/span>)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
Zooming in at x=1:<\/span><\/div>
plot(stiffSol.Time,stiffSol.Solution,<\/span>'-o'<\/span>)<\/span><\/span><\/div><\/div>
xlim([9620.2 9748.2])<\/span><\/span><\/div><\/div>
ylim([0.99464 1.00232])<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
Much better behaved! <\/span><\/div>
Being able to tell MATLAB that the problem is <\/span>\"stiff\"<\/span> makes life a little easier than before but I was a disappointed that <\/span>\"auto\" <\/span>didn't realise that my problem was stiff and choose a relevant solver for me. <\/span><\/div>
I reached out to development to ask why it failed me. They told me that the automatic solver selection is a heuristic that operates without peering inside your equations. It reacts to the data you supply to define the problem, e.g. if you've supplied a Jacobian, what the relative tolerances are, etc. There's no attempt to diagnose stiffness at all. That's why there are <\/span>\"auto\"<\/span>, <\/span>\"stiff\"<\/span>, and <\/span>\"nonstiff\"<\/span> automatic modes.<\/span><\/div>
Of course, this may change in future releases. Now this feature exists, it will be possible to improve on it and its already pretty useful! <\/span><\/div>

More flexible in time<\/span><\/h4>
Consider the following classic system of first-order ODEs that describes simple harmonic motion<\/span><\/div>
$ \\frac{dy_1}{dt} = y_2 \\ \\frac{dy_2}{dt} = -y_1 $<\/span><\/div>
subject to the initial conditions <\/span>y<\/mi><\/mrow>1<\/mn><\/mrow><\/msub>(<\/mo>0<\/mn><\/mrow>)<\/mo><\/mrow>=<\/mo>1<\/mn><\/mrow><\/mi><\/mrow><\/msub><\/mrow><\/math><\/span>and <\/span>y<\/mi><\/mrow>2<\/mn><\/mrow><\/msub>(<\/mo>0<\/mn><\/mrow>)<\/mo><\/mrow>=<\/mo>0<\/mn><\/mrow><\/math><\/span><\/div>
With the classic suite of solvers, you might solve this as follows<\/span><\/div>
y0 = [1; 0]; <\/span>% initial conditions: y1(0)=1, y2(0)=0<\/span><\/span><\/div><\/div>
tspan = [0 10]; <\/span>% time interval from 0 to 10 for the solution<\/span><\/span><\/div><\/div>
[t,y] = ode45(@(t,y)[y(2);-y(1)], tspan,y0);<\/span><\/span><\/div><\/div>
plot(t,y,<\/span>\"-o\"<\/span>)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
So far so simple. Imagine now that I want to go backwards in time from my initial condition as well as forwards. That is, I want my time span to be<\/span><\/div>
tspan2 = [-10 10]; <\/span>% time interval from 0 to 10 for the solution<\/span><\/span><\/div><\/div><\/div>
This is a little complicated using <\/span>ode45<\/span> since it assumes that the initial condition you set using <\/span>y0<\/span> corresponds to the beginning of the span. That is, if I do<\/span><\/div>
[t,y] = ode45(@(t,y)[y(2);-y(1)], tspan2,y0);<\/span><\/span><\/div><\/div>
plot(t,y,<\/span>\"-o\"<\/span>)<\/span><\/span><\/div><\/div>
title(<\/span>\"Not the boundary condition I wanted!\"<\/span>)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
it applies the boundary condition <\/span>y<\/mi><\/mrow>1<\/mn><\/mrow><\/msub>(<\/mo>-<\/mo>10<\/mn><\/mrow>)<\/mo><\/mrow>=<\/mo>1<\/mn><\/mrow><\/mi><\/mrow><\/msub><\/mrow><\/math><\/span>and <\/span>y<\/mi><\/mrow>2<\/mn><\/mrow><\/msub>(<\/mo>-<\/mo>10<\/mn><\/mrow>)<\/mo><\/mrow>=<\/mo>0<\/mn><\/mrow><\/math><\/span> which is not the problem I wanted it to solve! Instead I have to make two calls to <\/span>ode45 -- <\/span>one that goes forwards in time and the other that goes backwards.<\/span><\/div>
y0 = [1; 0]; <\/span>% initial conditions: y1(0)=1, y2(0)=0<\/span><\/span><\/div><\/div>
[tf,yf] = ode45(@(t,y)[y(2);-y(1)], [0 10],y0);<\/span><\/span><\/div><\/div>
[tb,yb] = ode45(@(t,y)[y(2);-y(1)], [0 -10],y0);<\/span><\/span><\/div><\/div>
plot([flipud(tb);tf],[flipud(yb);yf],<\/span>\"-o\"<\/span>)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
We got there but it took a couple of backflips. With the new interface, I don't need to worry about this as its much more flexible. <\/span><\/div>
shm = ode(ODEFcn=@(t,y)[y(2);-y(1)],InitialTime=0,InitialValue=[1;0],Solver=<\/span>\"ode45\"<\/span>);<\/span><\/span><\/div><\/div>
shmSol = solve(shm,-10,10);<\/span><\/span><\/div><\/div>
plot(shmSol.Time,shmSol.Solution,<\/span>'-o'<\/span>);<\/span><\/span><\/div><\/div>
title(<\/span>\"The same result but we got there more easily\"<\/span>)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
Unlike the traditional suite, the solution span we ask for doesn't even need to include our initial conditions. <\/span><\/div>
% Request a solution span that doesn't include the initial conditions at t=0<\/span><\/span><\/div><\/div>
shmSol2 = solve(shm,2*pi,4*pi)<\/span><\/span><\/div>
shmSol2 = <\/span><\/div>
ODEResults with properties:\n\n Time: [6.2832 6.4381 6.6484 6.8587 7.0690 7.2792 7.4806 7.6820 7.8833 8.0847 8.2869 8.4891 8.6913 8.8935 9.0955 9.2975 9.4995 9.7015 9.8989 10.0962 10.2935 10.4908 10.6933 10.8959 11.0984 11.3009 11.4952 11.6896 11.8839 12.0783 \u2026 ] (1\u00d734 double)\n Solution: [2\u00d734 double]\n<\/div><\/div><\/div><\/div><\/div><\/div>
The closeness of the solution points around t = 0 is the solver starting the two integrations with a small initial step size. If we would rather choose the output points ourselves, we can just specify them.<\/span><\/div>
shmSol3 = solve(shm,linspace(-pi,pi));<\/span><\/span><\/div><\/div>
plot(shmSol3.Time,shmSol3.Solution,<\/span>'-o'<\/span>);<\/span><\/span><\/div><\/div>
xlim([-pi,pi]);<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
We can get all of these results using the traditional suite, of course, it's just that its easier and more intuitive now. <\/span><\/div>

Event detection<\/span><\/h4>
Event detection has always been part of the MATLAB ODE suite and so, of course, this is also possible with the new framework. The <\/span>documentation for ode<\/span><\/a> demonstrates how to solve the classic bouncing ball problem using event detection. Here, I'll demonstrate how to use it to find maxima and minima of the solution to an ODE problem. The <\/span>original problem<\/span><\/a> comes from a 2011 blog post by John Kitchin.<\/span><\/div>
First, let's solve the ODE <\/span>y<\/mi><\/mrow>\u2032<\/mo><\/mrow><\/msup>=<\/mo>e<\/mi><\/mrow>-<\/mo>0<\/mn>.<\/mo>05<\/mn>t<\/mi><\/mrow><\/msup>sin<\/mi>(<\/mo>t<\/mi><\/mrow>)<\/mo><\/mrow><\/mrow><\/math><\/span><\/div>
myODE = ode(ODEFcn=@(t,y) exp(-0.05*t)*sin(t),InitialTime=0,InitialValue=0); <\/span><\/span><\/div><\/div>
sol = solve(myODE,0,20,Refine=10); <\/span><\/span><\/div><\/div>
plot(sol.Time,sol.Solution)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
We want to find the maxima and minima of this solution. We know from calculus that these occur when our original equation <\/span>y<\/mi><\/mrow>\u2032<\/mo><\/mrow><\/msup>=<\/mo>e<\/mi><\/mrow>-<\/mo>0<\/mn>.<\/mo>05<\/mn>t<\/mi><\/mrow><\/msup>sin<\/mi>(<\/mo>t<\/mi><\/mrow>)<\/mo><\/mrow><\/mrow><\/math><\/span>is equal to zero. <\/span><\/div>
odeEvent<\/span><\/a> objects define events that the solver detects while solving an ordinary differential equation. An event occurs when one of the event functions you specify crosses zero.<\/span><\/div>
So, we create an <\/span>odeEvent <\/span>object with an event function equal to our original ODE.<\/span><\/div>
E = odeEvent(EventFcn = myODE.ODEFcn)<\/span><\/span><\/div>
E = <\/span><\/div>
odeEvent with properties:\n\n EventFcn: @(t,y)exp(-0.05*t)*sin(t)\n Direction: both\n Response: proceed\n CallbackFcn: []\n<\/div><\/div><\/div><\/div><\/div><\/div>
This time, let's solve the ODE with this event definition<\/span><\/div>
myODE = ode(ODEFcn=@(t,y) exp(-0.05*t)*sin(t),InitialTime=0,InitialValue=0,EventDefinition=E); <\/span><\/span><\/div><\/div>
sol = solve(myODE,0,20,Refine=10)<\/span><\/span><\/div>
sol = <\/span><\/div>
ODEResults with properties:\n\n Time: [0 0.2000 0.4000 0.6000 0.8000 1 1.2000 1.4000 1.6000 1.8000 2 2.2000 2.4000 2.6000 2.8000 3 3.2000 3.4000 3.6000 3.8000 4 4.2000 4.4000 4.6000 4.8000 5 5.2000 5.4000 5.6000 5.8000 6 6.2000 6.4000 6.6000 6.8000 7 \u2026 ] (1\u00d7101 double)\n Solution: [0 0.0205 0.0797 0.1733 0.2969 0.4451 0.6122 0.7917 0.9766 1.1594 1.3319 1.4872 1.6201 1.7256 1.8002 1.8417 1.8489 1.8223 1.7635 1.6754 1.5622 1.4285 1.2800 1.1238 0.9664 0.8139 0.6717 0.5448 0.4377 0.3544 0.2983 \u2026 ] (1\u00d7101 double)\n EventTime: [2.8422e-14 3.1416 6.2832 9.4248 12.5664 15.7080 18.8496]\n EventSolution: [4.2572e-28 1.8503 0.2696 1.6201 0.4659 1.4521 0.6088]\n EventIndex: [1 1 1 1 1 1 1]\n<\/div><\/div><\/div><\/div><\/div><\/div>
The solution includes all of the times the event occurs. Let's plot them<\/span><\/div>
plot(sol.Time,sol.Solution,<\/span>\"-\"<\/span>)<\/span><\/span><\/div><\/div>
hold <\/span>on<\/span><\/span><\/div><\/div>
plot(sol.EventTime,sol.EventSolution,<\/span>\"ro\"<\/span>)<\/span><\/span><\/div><\/div>
hold <\/span>off<\/span><\/span><\/div><\/div>
xlabel(<\/span>\"Time\"<\/span>)<\/span><\/span><\/div><\/div>
ylabel(<\/span>\"y\"<\/span>)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
Exactly what we were looking for.<\/span><\/div>

Obtaining the ODE solution as a function<\/span><\/h4>
Until now, we have been working with discrete representations of the solutions to our ODE problems and this can present limitations on further analysis. For example, what if we wanted to integrate the solution? What we need is a function that represents the solution and we can get one using <\/span>solutionFcn<\/span><\/a>, passing in the <\/span>ode<\/span> object we are aiming to solve. <\/span><\/div>
myODE = ode(ODEFcn=@(t,y) exp(-0.05*t)*sin(t),InitialTime=0,InitialValue=0); <\/span>% Define the ODE problem<\/span><\/span><\/div><\/div>
solFcn = solutionFcn(myODE,0,20); <\/span>% Use it to create a solutionFcn <\/span><\/span><\/div><\/div><\/div>
I can now evaluate the solution at any point I like<\/span><\/div>
solFcn(15)<\/span><\/span><\/div>
ans = 1.3404<\/div><\/div><\/div><\/div>
or pass it to MATLAB's integral function<\/span><\/div>
integral(solFcn,0,20)<\/span><\/span><\/div>
ans = 19.5314<\/div><\/div><\/div><\/div>
Again, you could have done something similar with the old framework using <\/span>deval<\/span> and friends but I much prefer the new way! <\/span><\/div>

Mass Matrices, Jacobians, Parameters and everything else I haven't covered here<\/span><\/h4>
This is just an introduction to the new framework and I encourage you to <\/span>read the documentation<\/span><\/a> to explore everything else it offers. Let me know what you think in the comments section. <\/span><\/div>\n<\/div>