I am looking for a numerical method to help me solve the following system of two equations for $\theta(\xi,t)$ and $\lambda(\xi,t)$:
\begin{equation} \ddot{\theta}=\theta''+\lambda \\ \dot{\theta}^2\cos\theta+\ddot{\theta}\sin\theta=\lambda''+\lambda'\theta' \end{equation}
where $\dot{\theta} = \frac{\partial \theta}{\partial t}$ and $\theta'=\frac{\partial \theta}{\partial \xi}$.
I have the initial conditions $\theta(\xi,0)$ and $\dot{\theta}(\xi,0)$. I also have the boundary conditions $\theta(0,t)$, $\theta(\ell,t)$, $\lambda(0,t)$, and $\lambda(\ell,t)$.
I have tried to set this up for the builtin function pdepe to input into MATLAB, but I have not been able to put the equations into parabolic-elliptic form.
Thank you.
EDIT: The system stated above is a model system for me to try and get a hold of the problem. The equations I am actually trying to solve are a little different. They come from balances of linear and angular momentum for a rod which can bend. The rod has a torque, $T$, applied at its left end at time zero. I am trying to obtain the slope of the rod, $\theta$, and the internal axial and shear forces in the rod, $\lambda$ and $\mu$, for all values of space and time. The actual system is as follows:
\begin{equation} \ddot{\theta}=\theta''+\mu \\ \left( -\ddot{\theta}\sin\theta - \dot{\theta}^2\cos\theta \right)\cos \theta + \left( \ddot{\theta} \cos \theta - \dot{\theta}^2 \sin \theta \right) \sin \theta = \lambda'' - \lambda \theta'^2 - 2 \mu' \theta' -\mu \theta'' \\ \left( \ddot{\theta} \sin \theta + \dot{\theta}^2 \cos\theta \right) \sin \theta + \left( \ddot{\theta} \cos{\theta} -\dot{\theta}^2 \sin \theta \right) \cos \theta = \mu'' - \mu \theta'^2 + 2\lambda'\theta'+\lambda\theta''. \end{equation}
I have left out some of the quantities carrying physical dimension, such as the flexural rigidity and density, to focus on the mathematics of the problem.
The rod is initially straight and isn't moving. This gives me the IC's: $\theta \left( \xi,0 \right) = 0$ and $\dot{\theta} \left(\xi,0 \right)=0$.
There is a torque at the left end and no torque at the right end ($\xi = \ell$). This gives me the BC's: $\theta' \left(0^+,t\right) = T$ and $\theta' \left(\ell^-,t\right) = 0$.
The rod has no axial force at either the left or right end. This gives me the BC's: $\lambda \left( 0^+,t \right) = 0$ and $\lambda \left( \ell^-,t \right) = 0$.
The rod has no shear force at either the left or right end. This gives me the BC's: $\mu \left( 0^+,t \right) = 0$ and $\mu \left( \ell^-,t \right) = 0$.
One notes that there is no time derivative of $\lambda$. Substituting the first equation in the second one, we rewrite the system as \begin{aligned} \partial_t \theta &= \dot\theta \\ \partial_t\dot\theta &= \theta'' + \lambda \\ {\dot\theta}^2\cos\theta + (\theta'' + \lambda)\sin\theta &= \lambda'' + \lambda'\theta' \, . \end{aligned} Using finite differences on a regular grid $t_n = n\Delta t$, $\xi_i = i\Delta \xi$, one may write \begin{aligned} \theta_i^{n+1} &= \theta_i^n + \Delta t\,\dot\theta_i^n \\ \dot\theta_i^{n+1} &= \dot\theta_i^{n} + \Delta t\, \frac{\theta_{i+1}^n - 2 \theta_i^n + \theta_{i-1}^n}{{\Delta\xi}^2} + \Delta t \, \lambda_i^n \end{aligned} and \begin{aligned} &\left(\dot\theta_i^n\right)^2\cos\theta_i^n + \left(\frac{\theta_{i+1}^n - 2 \theta_i^n + \theta_{i-1}^n}{{\Delta\xi}^2} + \lambda_i^n\right)\sin\theta_i^n \\ =\; &\frac{\lambda_{i+1}^n - 2 \lambda_i^n + \lambda_{i-1}^n}{{\Delta\xi}^2} + \frac{\lambda_{i+1}^n - \lambda_{i-1}^n}{2\,{\Delta\xi}}\frac{\theta_{i+1}^n - \theta_{i-1}^n}{2\,{\Delta\xi}} . \end{aligned} The first two equations update $\theta$ and $\dot\theta$. The last equation provides the linear algebraic system \begin{aligned} &(1 + \tfrac{1}{4}(\theta_{i+1}^n - \theta_{i-1}^n))\lambda_{i+1}^n - (2 + {\Delta\xi}^2\sin\theta_i^n) \lambda_i^n + (1 - \tfrac{1}{4}(\theta_{i+1}^n - \theta_{i-1}^n))\lambda_{i-1}^n \\ =\; & {\Delta\xi}^2\left(\dot\theta_i^n\right)^2\cos\theta_i^n + (\theta_{i+1}^n - 2 \theta_i^n + \theta_{i-1}^n)\sin\theta_i^n \end{aligned} satisfied by $\lambda_1^n,\dots , \lambda_I^n$, which values are required in the time-update of $\dot\theta$. The initial conditions give $\theta_i^0 = \theta(\xi_i,0)$ and ${\dot\theta}_i^0 = \dot\theta(\xi_i,0)$. The boundary conditions give $\theta_0^n = \theta(0,t_n)$, $\theta_I^n = \theta(\ell,t_n)$, etc.