Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I can numercally solve a system of two ODEs for x and p:

dy1 = -y1 + y2
dy2 = -1 + y2

with the following Matlab code:

[t y]=rk4sys(20,0.05);
[t' y]

function [t,y] = rk4sys(n,h)
C = [0;-1];
M = [-1,1;0,1];
y(1, :) = [0, 0];
f =@(t,y)(C+M*y')';
t(1) = 0;
  for j = 1 : n
    k1 = f(t(j), y(j,:));
    k2 = f(t(j)+h/2, y(j,:)+h*k1/2 );
    k3 = f(t(j)+h/2, y(j,:)+h*k2/2 );
    k4 = f(t(j)+h, y(j,:)+h*k3);
    y(j+1,:) = y(j,:)+h*(k1+2*(k2+k3)+k4)/6;

The main program produces:

ans =

     0         0         0
0.0500   -0.0013   -0.0513
0.1000   -0.0050   -0.1052
0.1500   -0.0113   -0.1618
0.2000   -0.0201   -0.2214
0.2500   -0.0314   -0.2840
0.3000   -0.0453   -0.3499
0.3500   -0.0619   -0.4191
0.4000   -0.0811   -0.4918
0.4500   -0.1030   -0.5683
0.5000   -0.1276   -0.6487
0.5500   -0.1551   -0.7333
0.6000   -0.1855   -0.8221
0.6500   -0.2188   -0.9155
0.7000   -0.2552   -1.0138
0.7500   -0.2947   -1.1170
0.8000   -0.3374   -1.2255
0.8500   -0.3835   -1.3396
0.9000   -0.4331   -1.4596
0.9500   -0.4862   -1.5857
1.0000   -0.5431   -1.7183

Now, I have to change the element of (1,3) in the matrix above (that is 0(zero) for now) in order for the last number of the third column (that is -1.7183) to be equal to 0(zero).

In other words, the element (1,3), which is 0(zero) for now, is just my initial guess. I have to change this number to meet that y2 (at t = 21), that is -1.7183 for now, be 0(zero).

I can do this in Excel, but I have failed finding the way to do it in Matlab. Any help would be greatly appreciated.

share|cite|improve this question
up vote 4 down vote accepted

You are looking for a shooting method.

In other words, you have the following problem description: $$\mathbf{y}'(t) = f(\mathbf{y},t),$$ $$y_1(0) = 0,$$ $$y_2(\tau) = 0.$$

This is a boundary value problem, and you need to transform it into an initial value problem.

To do so, instead of solving the equations above, you are really solving $$\mathbf{y}'(t) = f(\mathbf{y},t),$$ $$y_1(0) = 0,$$ $$y_2(0) = b.$$

for some estimate of $b$. Your problem is to find $b$ such that the original condition is met.

One technique is known as shooting. You make an initial guess for $b$, compute the error at the right boundary for $y_2(\tau)$, and adjust $b$ based on your error. This is essentially a root-finding problem, and now you can use root-finding techniques to find the "root" $b$.

Essentially, the solution to an ODE initial value problem can be written as $y(t;y_0)$, or the solution $y(t)$ parameterized by the initial conditions $y_0$. This is because, broadly speaking, solutions to ODEs are not unique until you specify initial conditions. So the actual solution trajectory depends greatly on your choice of initial conditions $y_0$.

In your case, you have a second-order system, so you can compute $\mathbf{y}(t;\mathbf{y}_0)$, where $\mathbf{y}_0 = \begin{pmatrix} 0 \\ b\end{pmatrix}$. Since the trajectory that leads you to $y_2(\tau) = 0$ is unique, you can iterate over different values of $b$ to find the right result, and indeed, formulate the problem as solving $F(b) = 0$, where $F$ is your "function" that computes the solution to your ODE. In other words, treat the MATLAB code that you've written above as a function of the parameter $b$ that returns the value $y_2(\tau)$. Then, compute the root.

There are a number of ways to do this; the simplest method, although by far not the fastest, is through bisection. Newton's method may also work, but that is more difficult to set up. Broyden's method (the secant method) will also work.

EDIT: Below is a solution. If this is homework, you must attribute this code properly. This code is presented to show you how shooting works, not to solve the problem for you

The following code can be pasted in a single file called "odetest.m". This is not the highest-quality code; I banged it together to illustrate a point, and I tried to design it in such a way so as to be mathematical in how it was structured. Certainly it can be optimized. That's not the point.

First, we define an RK4 ODE solver. Then, we define a function $F(b)$ such that $F$ returns the value of $y_2(t=1)$. This value depends entirely on the parameter $b$; as such, we wish to find $b$ such that $F(b) = 0$.

We then use Broyden's method, which is a secant method, to compute the estimate for $b$ from an initial guess of both the value and the derivative. Finally, we call this, and then plot some things.

The resulting plot is here: Broyden's Method Shooting

Broyden's method is not the only way to solve this problem; however, it converges quite quickly and as you can see, it is easy to implement.

function odetest
    [i,lambda,lh] = broyden(@F,0,3);

    hold on
    colors = 'brgky';
    legend_string = {};
    for j=1:length(lh)
        legend_string{j} = ['j  = ' num2str(j) ', b = ' num2str(lh(j))];

    plot([0 1],[0 0],'k'); % Plot the x-axis


    answer = lambda

function x=F(b,ploton,color)
    % Extra arguments for plotting mode; when the function is called as a
    % function, just use a single argument
    ode = @(t,y)[-y(1)+y(2);-1+y(2)];
    y0 = [0; b];
    T = [0 1];
    h = 0.05;
    t = T(1):h:T(2);

    y = RK4(ode,.05,t,y0);
    if nargin > 1
    x = y(2,end);

function y=RK4(f,h,t,y0)
    n = length(t);
    y = zeros(length(y0),n);
    y(:,1) = y0;
    for i=1:n-1
        k1 = f(t(i),y(:,i));
        k2 = f(t(i)+h/2, y(:,i)+h/2*k1);
        k3 = f(t(i)+h/2, y(:,i)+h/2*k2);
        k4 = f(t(i)+h, y(:,i)+h*k3);
        y(:,i+1) = y(:,i)+h/6*(k1+2*k2+2*k3+k4);

function [i,lambda,lh] = broyden(func,IX,IJ)
% lambda = broyden(func,IX,IJ)
% IX -- initial guess for X
% IJ -- initial guess for J = dF/dX
% i -- number of steps needed to converge
% lambda -- solution
% lh -- lambda history, the computed lambda at each step
% Attempt to compute the solution to a system of nonlinear equations
% without explicitly requiring the computation of the derivatives of the
% function.
lh = [];
lambda = IX;
B = IJ;
    for i=1:30
        lh = [lh; lambda];
        f = func(lambda);
        s = B\(-f);
        x = s + lambda;
        y = func(x)-f;
        B = B + ((y-B*s)*s')/(s'*s);
        if norm(x-lambda)/norm(lambda) < 1e-10
            lambda = x;
        lambda = x;
share|cite|improve this answer
I am a young researcher in public policy. So, the word "homework" does not really apply to me. (It however reminds me of nice and young old days...) Thank you very much for your help. I will analyze your code and get back to you soon. Thanks again! – Bill TP Aug 15 '12 at 18:27
Dear Gorcenski, I believe you are a big figure here in math forum. Would you please take a look at the following problem?…. Someone already gave the Maple code. But I wonder if you could solve the problem in Matlab, not in a symbolic way. It is an old problem in city economics. Old but difficult to me. I heard that 20,30 years ago people used to solve it in Fortran -- OMG! I discussed this problem with one of my friends in engineering, but it was not satisfactory. Sad enough. Thank you very much in advance. – Bill TP Aug 15 '12 at 18:44
I am far from a big figure in any forum, I am simply someone trying to keep his mind fresh. I will take a look at that other question, but I cannot guarantee I will have the time to fully analyze it. – Emily Aug 15 '12 at 21:27
I appreciate it. I think the problem should interest you. It looks easy, but it seems not actually. I've known this problem for two years, and am still trying solving it via Matlab. – Bill TP Aug 15 '12 at 23:05

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.