function ode2joy()  
  % Define frequencies
       C=1;    D=2;    E=3;    F=4;    G=5;   A=6;
  f=[261.63  293.66  329.63  349.23  392.00 440.00];
  % Define tempo
  dt=1.5;     % time duration for a measure
  s=8000.;    % sample rate
  h=0.5*dt; q=0.25*dt; e=0.125*dt; t=0.375*dt;
  % Define score
  pitch   =[E E F G G F E D C C D E E D D E E F A A F E D C C D E D C C D D E C D E F E C D E F E D C D];
  duration=[q q q q q q q q q q q q t e h q q q q q q q q q q q q t e h q q q q q e e q q q e e q q q q];
  n= max(size(pitch));
  % Construct a vector with entire score
  score=[]; 
  for i=1:n
    omega=2*pi*f(pitch(i)); time=0:1./s:duration(i);
    y=sin(omega*time);
    score = [score y];    
  end
  % Play the score
  soundsc(score);
end
