It looks like time-dependent Schrodinger equation. So, you can introduce vector of solution Y = (y1(t), y2(t), y3(t)), matrix A(t) 3*3 in RHS which contains your ai, bi, ci functions ( dY/dt = AY), and then you can numerically solve it by approximating T-exp (time-ordered exponent). Divide time to very very small dt, and use Y(t+dt) = exp (A(t)dt)Y(t). Python has numpy and scipy, matrix exp is build function expm. Of cause, you may to know Y(t=0) to calculate Y(t).

Greetings to all who have contributed to the problem and I would like to say thanks to them for the time, enthusiasm, they have shown. 

For the time varying coefficients, the ODEs can be solved easily according to the functions that you have provided with but now the problem arises when the associated parameters of the ODEs not only time varying but also depend on other factors, say.

Example  ai(t,p) where p could be an array of parameters. Should I just modify the function definition?

I tried but I got some kind of error.

Could you be kind enough to help me out with that?

Thanking you in advance.

Best 

Thanks a lot for your help and the code, it was quite helpful I would say. But I would like to add the info that you might not have got it through my poor way of framing the question, apology for that.

So, I shall give it another try;

Just imagine that ai(t) was a time varying coefficient and now it depends not only on time but also say another set of parameters (like temperature, humidity, precipitation); so it should be like this ai(t, P, T, H);

 as an example the form of ai(t, P, T, H) could be like the following;

ai(t, P, T, H)= exp(T-15)*sin(P*t)+sin(H).

Please let me know if you could get it or not so that you can help me once again.

Appreciating your time and patience.

Best

You could use Modelica with OpenModelica and PyFMI for inclusion into Python code. Feel free to contact me for additional details.

Dear Suman,

It is possible to solve such a system of three ODEs in Python analytically, as well as being able to plot each solution. For example, if we wish to solve the following Predator-Prey system of ODEs:

dx/dt = - 0.5x

dy/dt = 0.5x - 0.4y

dz/dt = 0.4y

we would write the Python code (I use Spyder) as

import numpy as np

import matplotlib.pyplot as plt

from sympy import init_printing

init_printing()

integrate(x)

from sympy import Function, dsolve, Eq, Derivative, sin, cos, symbols

from sympy.abc import x, y, z, t

f = Function('f')

from IPython.display import display

eq = (Eq(Derivative(x(t),t),-0.5*x(t)), Eq(Derivative(y(t),t),0.5*x(t)-0.4*y(t)), Eq(Derivative(z(t),t), 0.4*y(t)))

display(dsolve(eq))

which yields the solutions:

x(t) = A Exp(-0.5t)

y(t) = 5A Exp(-0.5t) + B Exp(-0.4t)

and

z(t) = 4A Exp(-0.5t) - B Exp(-0.4t) + C

If you wanted to solve:

dx/dt = - xy2sin(t)

dy/dt = 5y sin(t)

then you need to replace (or add to) the last two lines of the above code with:

eq2 = (Eq(Derivative(x(t),t),x(t)*(y(t)**2)*sin(t)), Eq(Derivative(y(t),t),5*y(t)*sin(t)))

display(dsolve(eq2))

which returns the solutions:

x(t) = A Exp(B Exp(-10 cos(t)))

and

y(t) = C Exp(-5 cos(t))

I have just realised that you asked this question 2 years ago, so sorry if my answer comes too late.

Best,

Stephen Mason

Dear Suman,

Additionally, I thought you or anyone else viewing these postings would be interested in the following Python code, which should be easy for you to manipulate according to your own model.

This time, I have adopted the three coupled ODEs:

dx/dt = - y

dy/dt = 3z - x

dz/dt = 2x

and have written a code that returns, not only the general solutions for x, y and z, but also the respective particular solutions along with a plot of all three. The code is as follows:

import numpy as np

import matplotlib.pyplot as plt

from sympy import init_printing

init_printing()

integrate(t)

from sympy import Function, Indexed, Tuple, sqrt, dsolve, solve, Eq, Derivative, sin, cos, symbols

from sympy.abc import k, t

from sympy import solve, Poly, Eq, Function, exp

from sympy.abc import x, y, z, a, b, c

f = Function('f')

from sympy import Indexed, IndexedBase, Tuple, sqrt

from IPython.display import display

from sympy import *

from sympy.abc import *

from sympy.plotting import plot

init_printing()

t, C1, C2, C3 = symbols("t C1 C2 C3")

x, y, z = symbols("x y z", cls = Function, Function = True)

eq1 = Eq(Derivative(x(t),t), - y(t))

eq2 = Eq(Derivative(y(t),t), 3*z(t) - x(t))

eq3 = Eq(Derivative(z(t),t), 2*x(t))

eq4 = solve((eq1, eq2, eq3))

display(eq4)

print(eq4)

# Solutions with undetermined constants of integration, C1, C2 and C3:

soln = dsolve((eq1, eq2, eq3))

display(soln)

print(soln)

# Calculation of the above constants:

soln = dsolve((eq1, eq2, eq3), ics = {x: 3, y: 1, z: 2})

display(soln)

constants = solve((soln[0].subs(t,0).subs(x(0),3), soln[1].subs(t,0).subs(y(0),1), soln[2].subs(t,0).subs(z(0),2)),{C1,C2,C3})

# The following commands display both the general and particular

# solutions to the above code in conventional form:

xsoln = expand(soln[0].rhs.subs(constants))

display(xsoln)

print(xsoln)

ysoln = soln[1].rhs.subs(constants)

display(ysoln)

print(ysoln)

zsoln = soln[2].rhs.subs(constants)

display(zsoln)

print(zsoln)

# The following commands plot the particular solutions to the

# above code:

plot((xsoln, (t, 0, 3.1)), (ysoln, (t, 0, 3.1)), (zsoln, (t, 0, 3.1)))

The three particular solutions returned are:

x(t) = --10*sqrt(2)*exp(t)*sin(sqrt(2)*t)/11 + 25*exp(t)*cos(sqrt(2)*t)/11 + 8*exp(-2*t)/11

y(t) = 15*(2*sqrt(2)*sin(sqrt(2)*t) + cos(sqrt(2)*t))*exp(t)/11 + 5*sqrt(2)*(sin(sqrt(2)*t) - 2*sqrt(2)*cos(sqrt(2)*t))*exp(t)/11 + 16*exp(-2*t)/11

z(t) = 10*sqrt(2)*exp(t)*sin(sqrt(2)*t)/11 + 30*exp(t)*cos(sqrt(2)*t)/11 - 8*exp(-2*t)/11

Best regards

Stephen Mason