Implementing the $\theta$-method for solving an IVP

I want to implement the $\theta$-method to solve an IVP in MATLAB. The $\theta$-method is:

$y_{j+1} = y_j + h[\theta f(t_j, y_j) + (1 - \theta)f(t_{j+1}, y_{j+1})]$ for $\theta \in [0, 1]$.

I want to use it to solve an IVP of this form:

$y_1\prime = f_1(y_1, y_2) \;\;\;\;\; y_2\prime = f_2(y_1, y_2)$

$t \in [a, b]$

$y_1(a) = \alpha \;\;\;\;\; y_2(a) = \alpha \;\;\;\;\;\;$ (The initial conditions)

I'm not sure how I would implement this, since I'm used to implementing methods in which there is only one DE, such as $y\prime=f(t, y)$ (like Euler's method or the Runge-Kutta).

Here is what I did so far in MATLAB:

function thetaMeth(a, b, h, f, y)

j = 1;
theta = 0.5
for i = a:h:b
y(j+1) = y(j) + h*(theta*f(j, y(j)) + (1 - theta)*f(j+1, y(j+1)));
fprintf('t = %f; y(%d) = %f\n', i, j, y(j));
end
j = j + 1;
end


My parameter f would be for the $y_1\prime$ and $y_2\prime$ equations (it could possibly be taken in as a cell array of functions). I'm not sure what to do with the y parameter, because the solution equations $y_1$ and $y_2$ are not given (I'm not really sure if I should even have y as a parameter). This is an implicit method, and you need to have $y_{j+1}$ calculated beforehand somehow, because it appears in the LHS and RHS.

Anyone have any pointers as to how I can implement this? It doesn't have to be in MATLAB, I just need an algorithm that can be translated into code, since I'm not sure how to implement the algorithm myself.

-
You might also try posting this question on our sister site: scicomp.stackexchange.com – Potato Mar 11 '13 at 2:29
@Potato If I post it there, will I have to delete it here? – badjr Mar 11 '13 at 2:37
No. Just add a link to this one so people know that it has been cross-posted. – Potato Mar 11 '13 at 2:38
(1) If you have two differential equations then y should be a vector, i.e. you should have something like y(:,j+1) = y(:,j) + ... inside your loop where the first row corresponds to the values of $y_{1}$ and the second to the values of $y_{2}$. [Euler's method and Runge-Kutta can also be generalised in this way]. (2) Whenever you have an implicit method for solving a DE you will have to use fixed point iteration or a Newton-Rhapson type method at each step to find y(:,j+1). [Again, this should be familiar to you if you've implemented any implicit Runge-Kutta methods] – in_wolframAlpha_we_trust Mar 11 '13 at 9:07