Can Matlab solve differentail equations contain a solution from cubic equations?

Hi all, I'm fairly new to Matlab and only have a 2007 Matlab. I have the following equation that I'm trying to solve: And mf, xf,hf, ms, xs,and hs are all a function of time, t, such as xf=f(t)=xf(t),.

mf=xf*(0.2-hf-mf)^3/hf^3

nf=hf*(0.2-hf-mf)^3

ms=xs*(0.2-hs-ms)^3/hs^3

ns=hf*(0.2-hs-ms)^3

And I want to solve these 4 equations, which two of them are cubic equation. So I would get equations of mf and ms as funtions of xf(t), hf(t), xs(t), and hs(t) like:

mf=f(xf,hf)

nf=f(hf)

ms=f(xs,hs)

ns=f(hs)

And then I can solve differential equations using the solutions from above:

xf'[t]=mf-ms

xs'[t]=mf-ms

hf'[t]=hf-hs+3*xf[t]

hs'[t]=hf-hs+3*xf[t]

The initial conditions are:

xf[0]=0.001;

xs[0]=0;

hf[0]=0.001;

hs[0]=1;

[Roger Stafford]

Since your cubic equations can have up to three solutions, the initial conditions you have given may not be sufficient for solving the problem. You might have to specify which of the possible "branches" to follow.

In any case, cubic equations do have explicit solutions. For example, see

 http://en.wikipedia.org/wiki/Cubic_function

Selecting the appropriate solutions from these, you should be able to solve your differential equations using one of matlab's ode solvers.

I notice that nf and ns need not enter into your differential equations. They can be evaluated later.

Answer 1

[Walter Roberson]

Using dsolve() of the Symbolic toolbox, and eliminating contradictions,

nf(t) = -hf(t)*(hf(t)-1/5+mf(t))^3
xf(t) = -mf(t)*hf(t)^3/(hf(t)-1/5+mf(t))^3
ns(t) = -hf(t)*(-1/5+hs(t))^3

When you evaluate at t=0, everything is known in the xf equation except for mf, so this allows you to calculate mf(0), and substitute back in. mf(0) will have three solutions; one of those solutions leads to a contradiction, but the other two resolve completely (but if you calculate in the straight-forward way then you might seem to observe contradictions due to round-off)

There are possibly up to 4 other solutions that would need to be calculated a different way, as they give divide by 0 errors when t=0 is substituted in after the dsolve()