How do I use MATLAB to solve an Euler Backward problem for a Differential Equation? Is there a way to use MATLAB to solve an Euler Backward problem when the function I have is a differential equation? 
My original problem was to simulate an electromagnetic and an electrical field for a certain interval in R^2 and in this interval an electron would enter with an initial speed in x, y direction. I then plotted the path the electron would take inside the field, but now I'm trying to approximate the path using Euler Backward.
Our differential equation in MATLAB looks like this:
dy(1)=y(3);
dy(2)=y(4);
dy(3)=y(4)+1;
dy(4)=-y(3)

Where y=[x ; y ; u ; v] where u is the derivative of x and v is that of y.
 A: Given the vector-valued differential equation $\mathbf{y}'(t) = \mathbf{f}(\mathbf{y}(t))$, the Euler backward scheme is given by discretizing the time derivative and evaluating $\mathbf{f}$ at the value of $\mathbf{y}$ to be determined:
$$\frac{\mathbf{y}_{n+1} - \mathbf{y}_n}{h} = \mathbf{f}(\mathbf{y}_{n+1})$$
where $h$ is the time step and $t = nh$, so you have
$$\mathbf{y}_{n+1} - \mathbf{y}_n - h\mathbf{f}(\mathbf{y}_{n+1}) = 0$$
which is an implicit equation for $\mathbf{y}_{n+1}$. There are two interesting facts about the case you're trying to solve:
(1) the equation is basically linear, i.e. $\mathbf{f}(\mathbf{y}) = \mathbf{b} + A\mathbf{y}$ where $\mathbf{b}=(0,0,1,0)$ and $A$ is a matrix.
(2) the differential equation for $y_3$ and $y_4$ has no dependence on $y_1$ and $y_2$, implying that you can solve for $y_3$ and $y_4$ separately, and then calculate $y_1$ and $y_2$.
If we consider that $\mathbf{f}$ is linear, then we get the equation
$$\mathbf{y}_{n+1} - \mathbf{y}_n - hA\mathbf{y}_{n+1} - h\mathbf{b} = 0$$
which can be solved explicitly for $\mathbf{y}_{n+1}$, giving
$$\mathbf{y}_{n+1} = (I - hA)^{-1}(\mathbf{y}_n + h\mathbf{b})$$
where $I$ is the identity matrix. This requires finding the inverse of a 4 x 4 matrix which, whilst it is entirely possible using a computer, isn't necessarily desirable for stability reasons. It would be better to just consider the last two components of the equation, the $y_3$ and $y_4$ components. After manually inverting the 2 x 2 matrix you get
$$y_3^{n+1} = \frac{y_3^n + h y_4^n + h}{1+h^2}$$
$$y_4^{n+1} = \frac{y_4^n - hy_3^n}{1+h^2}$$
which you should easily be able to program in Matlab, and then use the values for $y_3$ and $y_4$ to calculate $y_2$ and $y_1$.
Note that I've given you a very specific implementation of backward Euler which is tailored for the problem at hand. A more general backward Euler solver would have to use a root-finding technique such as Newton-Raphson iteration to calculate $\mathbf{y}_{n+1}$ at each step.
