## 35 Reputation

6 years, 124 days

## Thanks@Mariusz Iwaniuk...

@Mariusz Iwaniuk  thanks for your help. but i still have some issues with the code i develop.

NULL
> restart;
> dsys := {diff(x[1](y), y) = 1, diff(x[2](y), y) = x[3](y), diff(x[3](y), y) = -x[6](y)-H[a]*x[3](y)-V[0]*x[3](y)-B[r]*x[8](y), diff(x[4](y), y) = x[5](y), diff(x[5](y), y) = -H[a]*P[m]*x[3](y)-V[0]*P[m]*x[5](y), diff(x[6](y), y) = x[7](y), diff(x[7](y), y) = -N[t]*x[7](y)^2-V[0]*P[r]*x[7](y)-N[b]*x[7](y)*x[9](y), diff(x[8](y), y) = x[9](y), diff(x[9](y), y) = -V[0]*S[c]*x[9](y)+N[t]^2*x[7](y)^2/N[b]+N[t]*P[r]*x[7](y)/N[b]+N[t]*x[7](y)*x[9](y), x[1](0) = 0, x[2](0) = 0, x[3](0) = 0, x[4](0) = 0, x[5](0) = 0, x[6](0) = 0, x[7](0) = -1, x[8](0) = 1, x[9](0) = -.68};
/
| d                d                      d
< --- x[1](y) = 1, --- x[2](y) = x[3](y), --- x[3](y) = -x[6](y)
| dy               dy                     dy
\

- H[a] x[3](y) - V[0] x[3](y) - B[r] x[8](y),

d
--- x[4](y) = x[5](y),
dy

d
--- x[5](y) = -H[a] P[m] x[3](y) - V[0] P[m] x[5](y),
dy

d                      d                         2
--- x[6](y) = x[7](y), --- x[7](y) = -N[t] x[7](y)
dy                     dy

- V[0] P[r] x[7](y) - N[b] x[7](y) x[9](y),

d                      d
--- x[8](y) = x[9](y), --- x[9](y) = -V[0] S[c] x[9](y)
dy                     dy

2        2
N[t]  x[7](y)    N[t] P[r] x[7](y)
+ -------------- + ----------------- + N[t] x[7](y) x[9](y),
N[b]              N[b]

x[1](0) = 0, x[2](0) = 0, x[3](0) = 0, x[4](0) = 0,

x[5](0) = 0, x[6](0) = 0, x[7](0) = -1, x[8](0) = 1,

\
|
x[9](0) = -0.68 >
|
/
> sys1 := eval(dsys, {B[r] = 1, H[a] = 5, N[b] = .1, N[t] = .1, P[m] = .8, P[r] = 10, S[c] = 1, V[0] = 1});
/ d                d
{ --- x[1](y) = 1, --- x[2](y) = x[3](y),
\ dy               dy

d
--- x[3](y) = -x[6](y) - 6 x[3](y) - x[8](y),
dy

d
--- x[4](y) = x[5](y),
dy

d
--- x[5](y) = -4.0 x[3](y) - 0.8 x[5](y),
dy

d
--- x[6](y) = x[7](y),
dy

d                        2
--- x[7](y) = -0.1 x[7](y)  - 10 x[7](y) - 0.1 x[7](y) x[9](y),
dy

d                      d
--- x[8](y) = x[9](y), --- x[9](y) = -x[9](y)
dy                     dy

2
+ 0.1000000000 x[7](y)  + 10.00000000 x[7](y)

+ 0.1 x[7](y) x[9](y), x[1](0) = 0, x[2](0) = 0, x[3](0) = 0,

x[4](0) = 0, x[5](0) = 0, x[6](0) = 0, x[7](0) = -1,

\
x[8](0) = 1, x[9](0) = -0.68 }
/
> sol := dsolve(dsys, numeric, method = classical[rk4], stepsize = 0.1e-1);
Warning, The use of global variables in numerical ODE problems is deprecated, and will be removed in a future release. Use the 'parameters' argument instead (see ?dsolve,numeric,parameters)
proc(x_classical)  ...  end;
> plots:-odeplot(sol, [[y, x[8](y)]], y = 0 .. 2, color = ["Red"]);
Warning, could not obtain numerical solution at all points, plot may be incomplete

 1 2 3 4 Page 4 of 4
﻿