Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

MATLAB Matlab and ODE's

  1. Jun 6, 2010 #1
    I'm looking at the following example (http://www.cyclismo.org/tutorial/matlab/control.html):

    y'=x^2-y^2, y(0)=1, where they use the Euler method to approximate numerical solutions.

    This is the code:

    I don't understand the third and fourth line at the top. What is [tex]y = 0*x[/tex] and how did they get y(1) = 1...?

    Last edited: Jun 6, 2010
  2. jcsd
  3. Jun 6, 2010 #2


    User Avatar
    Staff Emeritus
    Science Advisor
    Gold Member

    You want to make sure that the size of the 'y' vector is the same as the size of the domain vector x so you can plot one against the other one.. So they start with y=0*x which should make y a vector of all zeroes that is the same length as x. Setting y(1)=1 means that the first entry of y will be 1. This is the same thing as saying that y, at the value of x that is the first one recorded in the domain vector, should be 1 (so all they've done is chosen an initial condition)
  4. Jun 19, 2010 #3
    Ok. I want to put [tex]y''(x) + \left(\frac{2}{x}\right) y'(x) + y^n = 0[/tex]

    into the following template:

    Code (Text):
    function[y,t] = euler(fun,y0,t0,T,h)
    % [y,t] = euler(fun,y0,t0,T,h) -
    % This function computes the solution to the IVP y'(t) = fun(y,t),
    % y(t0)= y0 , for a given function "fun(t,y)" using Euler's method.
    % The function is defined either via an m-file fun.m, or using the
    % "inline" command. y0 is the initial value, T is the maximum time, h
    % is the stepsize and t0=initial time.
    % The output is a vector containing the approximate solution y_euler.
    y(1) = y0;
    t(1) = t0;
    for i=1:T/h
    y(i+1) = y(i) + h*feval(fun,t(i),y(i));
    t(i+1) = t(i) + h;
    % End of m-file euler.m
    There are a few difficulties. First, I am working with a second order differential equation. Secondly, I am interested in a continuous range of solutions for [tex]n[/tex]. For certain values of n, I have a transformation method that simplifies my equation by eliminating the y'(x). But, the main problem I'm dealing with is coding the "n" part. Whenever I try, an error message in the command box comes up that says 'matrix must be squared' and I don't know what this means.

    I would appreciate any help.

Share this great discussion with others via Reddit, Google+, Twitter, or Facebook