I have a differential equation in which the equation contains a function. I don't know how to input and call the function to generate the graph I want. can you help me fix it

function fv=fgliom(T,y)
p1=0.0068;
p2=0.012;
p3=0.002;   
i1=4.7*10^-8;
i2=4.7*10^-5;
i3=4.7*10^-8;
c1=510;
c2=510;
a1=510;
a2=510;
a3=510;
s1=1.8*10^-2;
s2=1.8*10^-3;
s3=1.8*10^-3;
v=2;
j=50;
b=0.2;
u=0;
fv=zeros(5,1);
fv(1)= p1*y(1)*(1-(y(1)/c1))-s1*y(1)*(y(2)+y(3))-(i1*y(1)*y(5))/a1+y(1);
fv(2)= p2*y(2)*(1-(y(2)+y(3))/c2)-s2*y(1)*y(2)-u*F[(y(5)]*y(1)-(i2*y(1)*y(5))/a2+y(2);
    fv(3)= p3*y(3)*(1-(y(2)+y(3))/c2)-s3*y(1)*y(3)+u*F*(y(5))*y(1);
    fv(4)= v*y(1)'*F*(-y(1)'/c1)*y(4)-(i3*y(4)*y(5))/a3+y(4);
    fv(5)= j-b*y(5);
    
    function fv=fnaik(T0,y0)
p1=0.0068;
p2=0.012;
p3=0.002;   
i1=4.7*10^-8;
i2=4.7*10^-5;
i3=4.7*10^-8;
c1=510;
c2=510;
a1=510;
a2=510;
a3=510;
s1=1.8*10^-2;
s2=1.8*10^-3;
s3=1.8*10^-3;
v=2;
j=50;
b=0.2;
u=10^-3;
fv=zeros(5,1);
fv(1)= p1*y0(1)*(1-(y0(1)/c1))-s1*y0(1)*(y0(2)+y0(3))-(i1*y0(1)*y0(5))/a1+y0(1);
fv(2)= p2*y0(2)*(1-(y0(2)+y0(3))/c2)-s2*y0(1)*y0(2)-u*F*(y0(5))*y0(1)-(i2*y0(1)*y0(5))/a2+y0(2);
fv(3)= p3*y0(3)*(1-(y0(2)+y0(3))/c2)-s3*y0(1)*y0(3)+u*F*(y0(5))*y0(1);
fv(4)= v*y0(1)'*F*(-y0(1)'/c1)*y0(4)-(i3*y0(4)*y0(5))/a3+y0(4);
fv(5)= j-b*y0(5);

function fv=fturun(T1,y1)
p1=0.0068;
p2=0.012;
p3=0.002;   
i1=4.7*10^-8;
i2=4.7*10^-5;
i3=4.7*10^-8;
c1=510;
c2=510;
a1=510;
a2=510;
a3=510;
s1=1.8*10^-2;
s2=1.8*10^-3;
s3=1.8*10^-3;
v=2;
j=50;
b=0.2;
u=10^-2;
fv=zeros(5,1);
fv(1)= p1*y1(1)*(1-(y1(1)/c1))-s1*y1(1)*(y1(2)+y1(3))-(i1*y1(1)*y1(5))/a1+y1(1);
fv(2)= p2*y1(2)*(1-(y1(2)+y1(3))/c2)-s2*y1(1)*y1(2)-u*F*(y1(5))*y1(1)-(i2*y1(1)*y1(5))/a2+y1(2);
fv(3)= p3*y1(3)*(1-(y1(2)+y1(3))/c2)-s3*y1(1)*y1(3)+u*F*(y1(5))*y1(1);
fv(4)= v*y1(1)'*F*(-y1(1)'/c1)*y1(4)-(i3*y1(4)*y1(5))/a3+y1(4);
fv(5)= j-b*y1(5);

clc;clear all;format long;
[T y]=ode45('ftuber',0 200,[10000 0 0 2000 0 500 140]',10^-7);
[T0 y0]=ode45('fnaik',0 200,[10000 0 0 2000 0 500 140]',10^-7);
[T1 y1]=ode45('fturun',0 200,[10000 0 0 2000 0 500 140]',10^-7);

figure (1);
plot(T,y(:,1),'.-',T0,y0(:,1),'.-',T1,y1(:,1),'.-');
title('Grafik A')
legend('u=50','u=30','u=25')
xlabel('waktu (hari)');
ylabel('MA(t) (sel)');

figure (2);
plot(T,y(:,1),'.-',T0,y0(:,1),'.-',T1,y1(:,1),'.-');
title('Grafik B')
legend('u=50','u=30','u=25')
xlabel('waktu (hari)');
ylabel('MA(t) (sel)');

figure (3);
plot(T,y(:,1),'.-',T0,y0(:,1),'.-',T1,y1(:,1),'.-');
title('Grafik C')
legend('u=50','u=30','u=25')
xlabel('waktu (hari)');
ylabel('MA(t) (sel)');

figure (4);
plot(T,y(:,1),'.-',T0,y0(:,1),'.-',T1,y1(:,1),'.-');
title('Grafik D')
legend('u=50','u=30','u=25')
xlabel('waktu (hari)');
ylabel('MA(t) (sel)');

figure (5);
plot(T,y(:,1),'.-',T0,y0(:,1),'.-',T1,y1(:,1),'.-');
title('Grafik E')
legend('u=50','u=30','u=25')
xlabel('waktu (hari)');
ylabel('MA(t) (sel)');

here i want

[Walter Roberson]

You indicate that F(x) = x>0 (in terms of MATLAB logical expression processing.)

Is that the same F as where you have

fv(2)= p2*y1(2)*(1-(y1(2)+y1(3))/c2)-s2*y1(1)*y1(2)-u*F*(y1(5))*y1(1)-(i2*y1(1)*y1(5))/a2+y1(2);

so for example did you intend F*(y1(5)) to handle what is shown in your diagram as ? If so then remember that * is multiplication not function application, so you would have needed F(y1(5)) not F*(y1(5))

However...

F(x) = 0 if x <= 0, F(x) = 1 if x > 0, is a discontinuous function. Which is a problem for the mathematics that is used for all of the ode*() functions: the way that they determine which points are suitable ones to estimate the curve, depends upon the first and second derivatives of all of the functions being continuous.

To deal with this, in ode45() you need to use event functions to detect the places where F(q(t)) or F(-q(t)) would change between 0 and 1, and you have to terminate the calculation there. Then you restart ode45 from the same place as you left off. When you use the ode*() routines, it is not exactly the case that you cannot use if statements or logical comparisons such as F = @(x) x>0 but it is the case that for any one call to ode45() the same branch needs to be used each time.

[Latifah Hanum]

thank you sir. I need to use event functions to detect the places where F(q(t)) or F(-q(t)) would change between 0 and 1, but i can't construct the script , pls help me

[Walter Roberson]

function [value,isterminal,direction] = eventfun(T1, y1)
  value = [y1(5)>0, -y1(1)>0];
  isterminal = [true, true];
  direction = [0, 0];  %crossing both ways
end

Answer 1

[Sam Chak]

Hi Ms. @Latifah Hanum,

Equation 6 is a kind of ON/OFF switch that can be constructed using a scaled version of the signum function.

x = -1:0.001:1;
F = sign(x)/2 + 0.5;
plot(x, F, 'linewidth', 1.5)
ylim([-0.5 1.5]), grid on

If g = y(1), s = y(2), r = y(3), n = y(4), q = y(5), then for

fq = sign(y(5))/2 + 0.5;

For , you can probably create a dummy function for Eq. (7)

gdot  = % type out Eq.(7) here
fgdot = sign(-gdot)/2 + 0.5;

[Walter Roberson]

fq = y(5) > 0;

is easier, faster.

... But both versions suffer from the discontinuity problem.

[Sam Chak]

Thanks @Walter Roberson for the valuable suggestion. 👍

You have advised OP in the above comment on how to apply ode45() when the more efficient approach fq = y(5) > 0 is used.

[Latifah Hanum]

sir @Sam Chak, can you please edit my script, cause there is part I don't understand yet

[Latifah Hanum]

function fv=fgliom(T,y)
p1=0.0068;
p2=0.012;
p3=0.002;   
i1=4.7*10^-8;
i2=4.7*10^-5;
i3=4.7*10^-8;
c1=510;
c2=510;
a1=510;
a2=510;
a3=510;
s1=1.8*10^-2;
s2=1.8*10^-3;
s3=1.8*10^-3;
v=2;
j=50;
b=0.2;
u=0;
fq = y(5) > 0;
gdot=p1*y(1)*(1-(y(1)/c1))-s1*y(1)*(y(2)+y(3))-(i1*y(1)*y(5))/a1+y(1);
fgdot = sign(-gdot)/2 + 0.5;

fv=zeros(5,1);
fv(1)= p1*y(1)*(1-(y(1)/c1))-s1*y(1)*(y(2)+y(3))-(i1*y(1)*y(5))/a1+y(1);
fv(2)= p2*y(2)*(1-(y(2)+y(3))/c2)-s2*y(1)*y(2)-u*fq*y(1)-(i2*y(1)*y(5))/a2+y(2);
fv(3)= p3*y(3)*(1-(y(2)+y(3))/c2)-s3*y(1)*y(3)+u*fq*y(1);
fv(4)= v*gdot*-(fgdot/c1)*y(4)-(i3*y(4)*y(5))/a3+y(4);
fv(5)= j-b*y(5);

that is true ? @Sam Chak

[Walter Roberson]

Earlier you had

fv(4)= v*y1(1)'*F*(-y1(1)'/c1)*y1(4)-(i3*y1(4)*y1(5))/a3+y1(4);

Is the y1(1)' intended to indicate complex conjugate transpose, or is it intended to indicate "prime" (that is, derivative) ?

It is confusing that your equation 7 defines but your equation 10 uses -- is the different than ?

[Latifah Hanum]

that is derivative. dg(t)/dt same with g'(t) right?

fv(4)= v*y1(1)'*F*(-y1(1)'/c1)*y1(4)-(i3*y1(4)*y1(5))/a3+y1(4);

what about it is true ?

fv(4)= v*gdot*-(fgdot/c1)*y(4)-(i3*y(4)*y(5))/a3+y(4);

[Sam Chak]

Hi @Latifah Hanum, take a look here.

Would advise you to check ALL equations again, one by one, term by term:

[Latifah Hanum]

function fv=fgliom(T,y)

p1=0.0068;

p2=0.012;

p3=0.002;

i1=4.7*10^-8;

i2=4.7*10^-5;

i3=4.7*10^-8;

c1=510;

c2=510;

a1=510;

a2=510;

a3=510;

s1=1.8*10^-2;

s2=1.8*10^-3;

s3=1.8*10^-3;

v=2;

j=50;

b=0.2;

u=0;

fq = sign(y(5))/2 + 0.5;

gdot=p1*y(1)*(1-(y(1)/c1))-s1*y(1)*(y(2)+y(3))-(i1*y(1)*y(5))/a1+y(1);

fgdot = sign(-gdot)/2 + 0.5;

fv=zeros(5,1);

fv(1)= p1*y(1)*(1-(y(1)/c1))-s1*y(1)*(y(2)+y(3))-(i1*y(1)*y(5))/a1+y(1);

fv(2)= p2*y(2)*(1-(y(2)+y(3))/c2)-s2*y(1)*y(2)-u*fq*y(1)-(i2*y(1)*y(5))/a2+y(2);

fv(3)= p3*y(3)*(1-(y(2)+y(3))/c2)-s3*y(1)*y(3)+u*fq*y(1);

fv(4)= v*gdot*-(fgdot/c1)*y(4)-(i3*y(4)*y(5))/a3+y(4);

fv(5)= j-b*y(5);

error :

>> fgliommm

Not enough input arguments.

Error in fgliommm (line 20)

fq = sign(y(5))/2 + 0.5;

where I should put the command