Solving an iterative, implicit Euler method in MATLAB

Question

I'm trying to solve an iterative problem that includes an implicit (backwards) Euler method to find successive time values for a given function. The numerical problem is shown here:

$$ \begin{aligned} \frac{u_{r,s+1}-u_{r,s}}{\Delta t}&+u_{r,s}\frac{u_{r+1,s+1}-u_{r-1,s+1}+u_{r+1,s}-u_{r-1,s}}{4\Delta x}\\ &\quad+\frac{\eta^*_{r+1,s+1}-\eta^*_{r-1,s+1}+\eta_{r+1,s}-\eta_{r-1,s}}{4\Delta x}\\ &=\frac{1}{3}(\alpha x)^2\frac{u_{r+1,s+1}-2u_{r,s+1}+u_{r-1,s+1}-u_{r+1,s}+2u_{r,s}-u_{r-1,s}}{\Delta x^2 \Delta t}\\ &\quad+\alpha^2x\frac{u_{r+1,s+1}-u_{r-1,s+1}-u_{r+1,s}+u_{r-1,s}}{\Delta x\Delta t}, \end{aligned} $$

where $s$ is the time iteration and $r$ is the location iteration, $x$ is the true location ($x=r\,\Delta x$), and $\alpha$, $\Delta x$, and $\Delta t$ are constants. The $\eta$ values are already solved, so I need to solve for $u(r,s+1)$. However, the solution includes values for $u(r-1,s+1)$ and $u(r+1,s+1)$, which requires this to be solved as a system of equations over $r$ for each $s$. I have boundary conditions of $u$ for all $r$ when $s=1$, as well as u at the first and last $r$ values for all $s$ (i.e. $u(1,s)$ and $u(r_\max,s)$). I'm not sure what to do from here though to solve for the "interior" $u$ values at $s+1$.

For reference, the original, non-discretized equation is here:

$$ \frac{\partial\overline{u}}{\partial t}+\overline{u}\frac{\partial\overline{u}}{\partial x}+\frac{\partial\eta}{\partial x}=\frac{1}{3}(\alpha x)^2 \frac{\partial^3\overline{u}}{\partial x^2\partial t}+\alpha^2 x\frac{\partial^2 \overline{u}}{\partial x \partial t}. $$

Answer 1

Thanks to TroyHaskin, I've realized that this is in reality a linear problem, since the u(s) values are known from the initial conditions, and is therefore easy to solve in Matlab. Thanks to everyone else for help as well.

EDIT: To be more specific, the system can be solved linearly by separating the u(s+1) terms and their coefficients from everything else. The solution takes the form [Aw,Ap,Ae]u = Q, where u = [u(r-1,s+1),u(r,s+1),u(r+1,s+1)]^T. Because this is a tridiagonal matrix, it can be solved with minimum storage (note that I have changed location r to i and time s to n in my code, just for personal consistency):

    A_E = zeros(N,1);
    A_W = zeros(N,1);
    A_P = zeros(N,1);
    Q = zeros(N+1,1);
    Q(1) = u(1,n+1);
    A_P(1) = 1;

    for i = 2:N
        A_W(i) = (-u(i,n)/(4*dx))-(((alpha*x(i))^2)/(3*(dx^2)*dt))+((alpha^2)*x(i)/(dx*dt));
        A_P(i) = (1/dt)+((2*((alpha*x(i))^2))/(3*(dx^2)*dt));
        A_E(i) = (u(i,n)/(4*dx)) - (((alpha*x(i))^2)/(3*(dx^2)*dt)) - ((alpha^2)*x(i)/(dx*dt));
        Q(i) = u(i,n)/dt - (u(i,n)/(4*dx))*(u(i+1,n)-u(i-1,n)) ...
            - (eta(i+1,n+1)-eta(i-1,n+1)+eta(i+1,n)-eta(i-1,n))/(4*dx) ...
            + (((alpha*x(i))^2)/(3*(dx^2)*dt))*(-u(i+1,n)+(2*u(i,n))-u(i-1,n)) ...
            + ((alpha^2)*x(i)/(dx*dt))*(-u(i+1,n)+u(i-1,n));
    end

        %Gaussian Elimination
        for i = 2:N
            A_P(i) = A_P(i) - (A_W(i)*A_E(i-1)/A_P(i-1));
            Q(i) = Q(i) - (A_W(i)*Q(i-1)/A_P(i-1));
        end

        %Back Substitution
        for i = N:-1:2
            u(i,n+1) = (Q(i) - (A_E(i)*u(i+1,n+1)))/A_P(i);
        end