Is there a way to plot a function f with equations (from a ODE system) inside of it?

Question

I am trying to solve this assignment for my Differential Equations class and I am encountering some problems.

Given this system of ODEs:

dy(1) = y(2) − y(3);
dy(2) = −3*y(2) − 5*sin(ω*y(1));
dy(3) = y(1)*y(2);

with omega=3. Initial values: y(1)=0; y(2)=4; y(3)=1;

Plot f between t=0 and t=20:

f(t) = |2*cos(y1(t)) + y2(t)|

I have first tried simulating the equation system with ODE45 with this code:

On editor:

function dy=modelo10(t,y)
    global omega

    dy = zeros(3,1);
    dy(1)= y(2) - y(3);
    dy(2)= -3*y(2) - 5*sin(omega*y(1));
    dy(3)= y(1)*y(2);

end

On command tab:

>>global omega;
>>omega = 3;
>>[T,Y] = ode45(@modelo10, [0,20], [0,4,1]); 
%% I assign the function to variable f
>>f = @(t) 2.*cos(y(:,1)) + y(:,2);
%% And finally plot it with the values t
>>fplot(f, [0,20]);

I get a blank graph with the following error:

Warning: Function behaves unexpectedly on array inputs. To improve 
performance, properly
vectorize your function to return an output with the same size and 
shape as the input
arguments. 
  In matlab.graphics.function.FunctionLine>getFunction
  In matlab.graphics.function.FunctionLine/updateFunction
  In matlab.graphics.function.FunctionLine/set.Function_I
  In matlab.graphics.function.FunctionLine/set.Function
  In matlab.graphics.function.FunctionLine
  In fplot>singleFplot (line 241)
  In fplot>@(f)singleFplot(cax,{f},limits,extraOpts,args) (line 196)
  In fplot>vectorizeFplot (line 196)
  In fplot (line 166) 
Warning: Error updating FunctionLine.

The following error was reported evaluating the function in 
FunctionLine update:
Non-scalar in Uniform output, at index 1, output 1.
Set 'UniformOutput' to false.

Warning: Error updating FunctionLine.

The following error was reported evaluating the function in 
FunctionLine update:
Non-scalar in Uniform output, at index 1, output 1.
Set 'UniformOutput' to false

So my question is what are the steps (or commands) to plot f?

I am using MATLAB R2019a.

Answer 1

global omega;
omega = 3;

[T,Y] = ode45(@modelo10, [0,20], [0,4,1]);

• From here Y is a function of T --> No need to define new function Y = Y(T) for T = [0,20], this is no longer variable but real numbers Y = [Y1(0),..., Y1(20); Y2(0),..., Y2(20); Y3(0),..., Y3(20)]

• I assign the function to variable f

  f = @(t) 2.*cos(y(:,1)) + y(:,2) --> f = abs( 2.*cos(Y(:,1)) + Y(:,2));
  

  lowercase y and UPPERCASE Y are two different variables. Absolute value function in matlab is abs

• And finally plot it with the values t

  Use plot instead of fplot
  

  Actually you don't need to specify input value for function f as it has been already included when calculating ode45 whereas fplot takes a key argument, namely, function expression

  (function handle --> function defined as follow f = @(x)_____). 
  

  The input could be additionally added as second argument.

  It's optional; by default fplot evaluates function in the range [0 5]

• Here f is real number not function expression, it has been already evaluated

      plot(T, f);
  

Function f Plot for T ranging from 0 to 20

Answer 2

Your first problem is in defining f:

[T,Y] = ode45(@modelo10, [0,20], [0,4,1]); 
f = @(t) 2.*cos(y(:,1)) + y(:,2);

Note that the first line defines a variable Y, but in the second one you use y. These are different variables. This likely works for you because you have a variable y defined in your workspace, but it is doing the wrong thing.

Next, your f is a function f(t), but it is defined not using the input t. It always outputs the same thing. fplot expects a function that outputs an array of the same size as its input, this is the error message you are getting.

However, you don't need to define a function here, you can directly compute all values for f and plot them using plot:

[t,y] = ode45(@modelo10, [0,20], [0,4,1]);
f = 2.*cos(y(:,1)) + y(:,2);
plot(t,f)

---

Additionally, I want to recommend that you don't use global. You don't need it here, you can define omega as an input argument to modelo10:

function dy = modelo10(t,y,omega)
   dy = zeros(3,1);
   dy(1) = y(2) - y(3);
   dy(2) = -3*y(2) - 5*sin(omega*y(1));
   dy(3) = y(1)*y(2);
end

and then call ode45 as follows:

[t,y] = ode45(@(t,y)modelo10(t,y,omega), [0,20], [0,4,1]);

Here, @(t,y)modelo10(t,y,omega) defines an anonymous function that carries the value of omega with it. This anonymous function has two input arguments (t and y), as required by ode45.

Finally, you can simplify your code by defining modelo10 on a single line:

function dy = modelo10(t,y,omega)
   dy = [ y(2) - y(3); -3*y(2) - 5*sin(omega*y(1)); y(1)*y(2) ];
end

You can now do:

modelo10 = @(t,y) [ y(2) - y(3); -3*y(2) - 5*sin(omega*y(1)); y(1)*y(2) ];
[t,y] = ode45(modelo10, [0,20], [0,4,1]);