# Is it possible to solve or approximate this second order nonlinear system of differential equations.?

Given initial values $d[0]$ and $k[0]$, I would like to solve for the initial rate of change, $\dot d[0]$, and compare this value against some data.

I have the following profit function, which I would like to maximize:

$$\pi=\int^T_0 e^{-pt}(-50 (d - .1)^2 (d - .95) (d - .7) - 4 (.3 - k)^2 + 2 d k - \dot d^2-\dot k^2)$$

Where d and k are functions of time, t. I get the following Euler equations:

$$4020 d - 11100 d^2 + 8000 d^3 + 80 p \dot d = 299 + 80 k + 80 \ddot d$$ $$6 + 5 d + 5 \ddot k =20 k + 5 p \dot k$$

And the following boundary values: $$d[0]=\alpha, k[0]=\beta$$ $$-2 e^{-pT} \dot d[T] = 0,-2 e^{-pT} \dot k[T] = 0$$ Ideally, I would like an analytic solution in terms of $\alpha, \beta, p,$ and $T$ but specifying them from the start would also be ok.

I have tried solving the system in Mathematica with no luck. DSolve just returns the input, NDSolve complains about stiffness and singularity. I don't know how to deal with this. I successfully used a Taylor expansion to do approximate $\dot d[0]$, but given that there are supposed to be two stationary points I'm not sure how I can incorporate both so that I can test my predictions against my fairly large data set. How could I non-arbitrarily determine which stationary point to approximate the behavior around given $d[0]$ and $k[0]$? Also, I'm not sure how appropriate the approximation is, given that most of the values I'm interested in will be somewhat far from the stationary points.

Is there any way at all that I can solve or approximate a value for $\dot d[0]$, either analytically or numerically?

Any help would be greatly appreciated!

You could try what is called a shooting method. You can imagine values for $\dot d(0), \dot k(0)$, solve numerically and obtain $\dot d(T), \dot k (T)$. You can try to approximate the correct $\dot d (0), \dot k (0)$ so that it gets the right value.
To get a numerical code mesh an square $[-M,M] \times [-L,L]$. Then run a numerical code for your system of equations for each $(\dot d (0), \dot k(0)$ of your mesh. If you have a solution of your system with such $(\dot d (0), \dot k (0)) \in [-M,M] \times [-L,L]$ then you should get some value of $(\dot d(T), \dot k(T))$ close to the one you want. This comes as consequence of the continuous dependence of solutions with respect to the initial data.