Finds all odd pseudoprimes to base 2 which are <= m for a value of m input by the user

% Finds all odd pseudoprimes to base 2 which are <= m for a value of m input by the user.

m=input('Enter the largest number to be tested: ');
disp('Pseudoprimes to base 2 are: ')
for n=3:2:m % Want to test odd numbers n from 3 to m.
if pow(2,n-1,n)==1 % This tests whether 2^(n-1) leaves remainder
% 1 on division by n.
% If it doesn't then we go on to the next n.
% If it does then we test whether n is in fact
% prime by trial dividing n by odd numbers
% less than the square root of n. If none of
% these divides exactly into n then n is
% certainly prime, so should not be included
% in the list of pseudoprimes.

k=3; % The first odd number to trial divide into n
while k<=sqrt(n)
if rem(n,k)==0 % Tests whether k divides exactly into n
k=n+1; % This just gets out of the while loop,
% because n has been shown to be composite
% (multiple of k) so can be displayed
% as a pseudoprime.
else k=k+2; % otherwise go on to the next odd k to trial
% divide into n.
end;
end;
if k==n+1 disp(n); end;
end; % of if pow==1 loop
end; % of for n=3:2:m loop

Power algorithm for finding rem(a^n,m)

% This implements the power algorithm for finding rem(a^n,m) (the remainder on dividing a^n by m) reliably for large numbers (up to about 8 digits).
%
% Usage: Type pow(a,n,m) where a,n,m are already defined, or are explicit numbers as in pow(2,1000,25).


function x=pow(a,n,m)
b = a;
x = 1;
while n > 0
d=rem(n,2);
if d==1
x=rem(x*b,m);
end
b=rem(b * b,m);
n=(n-d)/2;
end

M-file to find primes < 5000

% M-file to find primes < 5000

r=11;
p=[2,3,5,7];
i=4; %i counts the number of primes so far
while r<5000
s=3;
while s*s<=r & rem(r,s)>0
s=s+2;
end;
if rem(r,s)>0
i=i+1;
p(i)=r;
end;
r=r+2;
end;

Cycloid

% Approximate times of descent for a cycloid from (0,c) to (d,0), corresponding to theta=0, theta=theta1 respectively and the value r for the radius of the rolling circle.
% Uses quad8 for theta=.01 to theta1 and estimates the integral from 0 to .01 separately.

g=0.1; % g is "gravity". 0.1 seems to produce reasonable values for the time

global r_gl c_gl

c_gl=input('Type c ');
d=input('Type d ');

% The values of r and theta you need here are calculated as
% in the text, using the M-file cycfun.m and the command fsolve.

r_gl=input('Type r ');
theta1=input('Type theta1 ');

time=quad8('slide4fn',.01,theta1);
% This uses Matlab's numerical integrator to find the
% time from theta=.01 to theta=theta1.

% Now we estimate the time from theta=0 to theta=.01

theta=.01

%******************************************************************
% These are the formulae for x, y and their derivatives
% x1, y1 with respect to theta.

x=r_gl*(theta-sin(theta));
y=c_gl-r_gl+r_gl*cos(theta);
x1=r_gl*(1-cos(theta));
y1=-r_gl*sin(theta);

%******************************************************************

initial=sqrt(1+(x1/y1)^2)*2*sqrt(c_gl-y)/sqrt(2*g)

time=time+initial;

time

Function snowboat

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%snowboat.m%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% function snowboat
%
% Solver for 1st order equations associated with snowboat problem
%
% Solution is for Vector [u v x y]' where u=x' , v=y'
% It prompts for:
% mean coefft. of kinetic friction mu
% coefft of air resistance kr
% starting coords [x0 y0]
% starting veloc. [u0 v0]
% time interval dt
% and requests permission to complete another time step.
% The derivatives specifying the motion are given in snbtfn.m
%
% Uses: snbtfn.m, snowsl.m, fsnow1.m
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function snowboat
% Initialization
hold off; clg;format bank;
disp('Equations of motion for snow boat')
disp(' ')
disp('When you press RETURN you will see a contour plot and then')
disp(' an empty plan view of the area')
disp('You will then be asked for the initial coordinates/velocities')
disp('Press RETURN to end any pauses.')
pause;
next = 'y';
dt = 0.5; trace = 0; tol = 0.05;
%
a=2000; nx=20; dx=a/nx; % square area 2km x 2km
[X,Y]=meshdom(0:dx:a,0:dx:a); % set up contour plot of the hill
Z=snowsl(X,Y);
contour(Z,20)
title('Countour plot of snow slopes')
pause;
clg;
%
xp = [0,2000]; yp = xp;

plot(xp,yp,'.r')
title('Plan view of snow slopes')
xlabel('x')
ylabel('y')
hold on
% pause

% Input parameters
mu = input(' Input coefft of friction : ');
kr = input(' Input coefft of air res. : ');
m = 150.0; % mass of snowboat + one passenger (kgs)
k= [mu,kr,m];
while next ~= 'n'
initpos = input('Initial x,y position, enter as [x0 y0]: ');
initvel = input('Initial velocity, enter as [u0 v0]: ');
dt=input('Suggest a time step (try 10 secs if you have no better suggestion): ');
disp('Initial z and vertical velocity component are');
z0=snowsl(initpos(1),initpos(2)); [fx fy]=fsnow1(initpos(1),initpos(2));
vz0=fx*initvel(1)+fy*initvel(2);
disp([z0 vz0])

disp('please wait...')
y0 = [initpos,initvel]'; t0 = 0; tf = dt; step = 'y'; mti=1;
while step == 'y'
[tout,yout] = ode23k('snbtfn', k, t0, tf, y0, tol, trace);
%
xp = yout(:,1); yp=yout(:,2); % these are output x,y, coords
% xp = yout(:,3); yp=yout(:,4); % these are output veloc. components
%
np = length(xp); mtf=mti+np-1;
plot(xp,yp,'r');
% pause

% plot this segment of run
disp(' t x(t) y(t) u(t) v(t)')
disp([tf xp(np) yp(np) yout(np,3) yout(np,4)])
disp(' z(t) vz(t)');
z=snowsl(xp(np),yp(np)); [fx fy]=fsnow1(xp(np),yp(np));
vz=fx*yout(np,3)+fy*yout(np,4);
disp([z vz])

% store the height and speed as a function of distance from starting point
yp1(mti:mtf)=snowsl(xp,yp); %height
yp2(mti:mtf)=10.0*sqrt(yout(:,3).^2+yout(:,4).^2); %speed
xp=xp-initpos(1); yp=yp-initpos(2);
xpp(mti:mtf)=sqrt(xp.*xp+yp.*yp); %horiz. dist.
mti=mtf+1;

step = input('Another time step? (y/n): ','s');
t0 = tf; tf = t0 + dt;
num = size(tout); y0 = yout(num(1),:)';
end
next = input('Do you want another simulation? (y/n): ','s');
end

% now plot the height and speed overall as promised
hold off
clg
disp('Press ENTER for a plot of the height Z and the speed V ...');
disp(' as a function of the horiz. dist. travelled ');
pause

xp = [0,1500]; yp = [400,800];
plot(xp,yp,'.r')
title('Height Z is White, Speed V (x10+400) is Green')
xlabel('d')
ylabel('Z,V')
hold on
plot(xpp,yp1,'w')
yp2=yp2+400; % add 400 to make visible on same plot
plot(xpp,yp2,'g')

Species simulates the behaviour of an ecological system over a time interval

function species
%
% species simulates the behaviour of an ecological system over a time interval 'dt'। It is based on the classic Lotka-Volterra model.

% It prompts for starting populations 'X0' and 'Y0,' and permission to complete another time step.
%
% species.m uses: ode23k.m, specfn.m
%

% Initialization
disp('Ecological system with two competing species.')
disp(' ')
disp('When you press RETURN you will see an empty plot of #Y vs #X.')
disp('You will be asked for the model parameters and an initial population.')
disp('Press RETURN to end any pauses.')
pause
more = 'y';
dt = 0.1; trace = 1; tol = 0.1;
xp = [0 10]; yp = [0 10];
k = [0 0 0 0];
clg; hold off
plot(xp,yp,'.')
title('Autonomous system of competing species X and Y')
xlabel('#X')
ylabel('#Y')
hold on
pause

% Input parameters
kx = input(' Model parameters for dx/dt eqns, enter as [a b c] : ');
ky = input(' Model parameters for dy/dt eqns, enter as [d e f] : ');
k(1:3)=kx; k(4:6)=ky;
while more ~= 'n'
pop0 = input('Initial populations (from 0 to 24), enter as [X0 Y0]: ');
dt=input('Suggest a time step (try 0.1 if you have no better suggestion): ');
disp('please wait...')
y0 = pop0'; t0 = 0; tf = dt; step = 'y';
while step == 'y'
[tout,yout] = ode23k('specfn', k, t0, tf, y0, tol, trace);
xp = yout(:,1); yp=yout(:,2);
plot(xp,yp,'r'); pause
step = input('Another time step? (y/n): ','s');
t0 = tf; tf = t0 + dt;
num = size(tout); y0 = yout(num(1),:)';
end
more = input('Do you want another simulation? (y/n): ','s');
end
hold off

Load in tide time and height data

%%%%%%%%%%%%%%%%%%%%%%tides.m%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% load in tide time and height data
%
%
% Usage:-
% tides
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
echo off
load times.mat
load heights.mat
%
disp('Tide times and heights read in OK')
disp('Times in tm(1:nt), heights in hm(1:nt)')
disp(' nt dt tmin tmax ')
format bank;
nt=length(tm); dt=tm(2)-tm(1);
disp([nt,dt,tm(1),tm(nt)])
format short;

To show, for the `tomato' data, the effect of finding a line fit

% M-file = toms.m

%

% To show, for the `tomato' data, the effect of finding a line fit (not sorting data), then aquadratic fit (not sorting data) and finally a quadratic fit (sorting data). Clearly for a quadratic fit, the sorting is essential! For the linear fit it is not essential because the line, though drawn as a zigzag, still looks like a line.



tomato

C=polyfit(x,y,1);

Y=polyval(C,x);

plot(x,y,'*w',x,Y,'w')

xlabel('This is the linear fit, without sorting')



pause



C=polyfit(x,y,2);

Y=polyval(C,x);

plot(x,y,'*w',x,Y,'w')

xlabel('This is the quadratic fit, without sorting!')



pause



[sx i]=sort(x);

sy=y(i);

C=polyfit(sx,sy,2);

Y=polyval(C,sx);

plot(sx,sy,'*w',sx,Y,'w')

xlabel('This is the quadratic fit, WITH sorting')

Produce and plot a nice topographical surface

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%topog.m
%
% Produce and plot a nice topographical surface
%
% Then allows one find a local minimum of the surface
%
% Usage:-
% topog
%
% then follow the prompts.
%
% uses: snowls.m, snowmn.m, snowmx.m
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
echo off
clg;
a=2000; nx=20; dx=a/nx;
[X,Y]=meshdom(0:dx:a,0:dx:a);

Z=snowsl(X,Y);

mesh(X,Y,Z)

hold on
title('Snow slopes')
next=input(' More? (y/n): ','s');
if next == 'y',
hold off;
clg
contour(Z,20)
X0=input(' Input guess for Minimum: [x,y]: ');
disp('please wait...');
format bank;
Xmin=fmins('snowmn',X0);
result=[Xmin,snowsl(Xmin(1),Xmin(2))];
disp(' Xmin Ymin Zmin');
disp(result);
X0=input(' Input guess for Maximum: [x,y]: ');
disp('please wait...');
Xmax=fmins('snowmx',X0);
result=[Xmax,snowsl(Xmax(1),Xmax(2))];
disp(' Xmax Ymax Zmax');
disp(result);
end

Draws the Taylor approximation to a sine curve of degree 2k-1

% Draws the Taylor approximation to a sine curve of degree 2k-1
k=input('Type the number of terms of the Taylor series for sine: ');
x=0:.05:2*pi;
z=sin(x);
plot(x,z,'g') %Actual sine curve is drawn in green
hold on
pause % This holds up execution until Enter is pressed
w=x;
y=x;
s=-1;
for j=1:k-1
w=w.*x.*x/(2*j*(2*j+1));
y=y+s*w;
s=-s;
end
plot(x,y,'r') % approximation is drawn in red
hold off

Draws the Taylor approximation to a sine curve or order n

% Draws the Taylor approximation to a sine curve or order n but here the scaling is left fixed after the actual sine curve is drawn

k=input('Type the number of terms of the Taylor series for sine: ');

x=0:.05:6.30;

z=sin(x);
plot(x,z,'g') %Actual sine curve is drawn in green

axis(axis) % Freeze the scaling
hold on

pause

y=x;
s=-1;
w=x;
for j=1:k-1
w=w.*x.*x/(2*j*(2*j+1));
y=y+s*w;
s=-s;
end;
plot(x,y,'r') %Approximation is drawn in red
hold off

vderpol solves the second order Van der Pol equation

function vderpol
%
% vderpol solves the second order Van der Pol equation
%
% w" -k*(1-w*w)*w' + w = 0
%
% as an autonomous first order system in v (=w') and w
% over a time interval 'dt'. It prompts for starting
% coordinates 'w0' and 'v0,' and requests
% permission to complete another time step.
% The derivatives specifying the equation are given in
% vdplfn.m
%
% vderpol.m uses: vdplfn.m, ode23k.m
%

% Initialization
disp('First order system for Van der Pol equation.')
disp(' ')
disp('When you press RETURN you will see an empty plot of v vs w.')
disp('You will be asked for the model parameter k, and initial coordinates.')
disp('Press RETURN to end any pauses.')
pause
more = 'y';
dt = 0.5; trace = 0; tol = 0.05;
xp = [-5 5]; yp = [-5 5];
k = 0;
clg; hold off
plot(xp,yp,'.')
title('1st order system for Van der Pol equation')
xlabel('w')
ylabel('v')
hold on
pause

% Input parameters
k = input(' Model parameter k : ');
while more ~= 'n'
init = input('Initial coordinates (from -4 to 4), enter as [w0 v0]: ');
dt=input('Suggest a time step (try 0.5 if you have no better suggestion): ');
disp('please wait...')
y0 = init'; t0 = 0; tf = dt; step = 'y'; mti=1;
while step == 'y'
%
[tout,yout] = ode23k('vdplfn', k, t0, tf, y0, tol, trace);
%
xp = yout(:,1); yp=yout(:,2);
nti = length(tout); mtf=mti+nti-1;
xpp(mti:mtf)=tout(1:nti);
ypp(mti:mtf)=xp(1:nti); mti=mtf+1;
plot(xp,yp,'r');
pause
disp(' t w(t) v(t)')
disp([tf xp(nti) yp(nti) ])
step = input('Another time step? (y/n): ','s');
t0 = tf; tf = t0 + dt;
num = size(tout); y0 = yout(num(1),:)';
end
more = input('Do you want another solution? (y/n): ','s');
end
hold off
clg
disp('Press ENTER for a plot of the last solution: w(t) vs t ');
pause
plot(xpp,ypp,'g');
title('Van der Pol equation')
xlabel('t')
ylabel('w')

Draws a zigzag based on lengths of 100

% Draws a zigzag based on lengths of 100 and l and angles theta1, theta2 (input in degrees).
% This is the most basic version where you have to provide the size of the graphics window and the number of steps.

l=input('Type the value of l ');
theta1=input('Type the angle theta1 in degrees ');
theta2=input('Type the angle theta2 in degrees ');
steps=input('Type the number of steps ');
xl=input('Type the lower bound of x ');
xu=input('Type the upper bound of x ');
yl=input('Type the lower bound of y ');

yu=yl+xu-xl; % This is to make the window square

theta1=theta1*pi/180; % Convert to radians!
theta2=theta2*pi/180;



x1=100; y1=0; x2=100+l; y2=0;
% (x2,y2) is the point reached at the end of step 1

x=[0 x1]; y=[0 y1]; plot(x,y,'g');
% This plots the line joining (0,0) to (100,0) in green

axis([xl,xu,yl,yu]);
axis('square');


hold on % So that all plots are made on the same picture

x=[x1 x2]; y=[y1 y2]; plot(x,y,'r');
% This plots the line joining (100,0) to (100+l,0) in red

for j=1:steps-1 % Now we do the zigzag, starting with step 2

x1=x2+100*cos(j*theta1);
y1=y2+100*sin(j*theta1);
% These calculate the end of the next 'zig'

x=[x2 x1]; y=[y2 y1]; plot(x,y,'g');
% This plots the next 'zig' in green

x2=x1+l*cos(j*theta2);
y2=y1+l*sin(j*theta2);
% These calculate the end of the next 'zag'

x=[x1 x2]; y=[y1 y2]; plot(x,y,'r');
% This plots the next 'zag' in red
end;

hold off

Loads a noisy signal similar to fftdemo.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%fourdat.
%
% Loads a noisy signal similar to fftdemo.m. Same as that used by
% fouran.m
%
% Usage:-
% fourdat
%
% and press return as required
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

echo off
% load the data
load noisyt.mat;
load noisys.mat;

disp('Press RETURN for a plot of some of the data')
pause
plot(t(1:50),s(1:50),'g')

load in signal data and then Fourier analyse

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% fouran.m
%
% load in signal data and then Fourier analyse
%
%
% Usage:-
% fouran
%
% press return as required
%
% 'fourdata' loads the same data with no Fourier analysis
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function fouran

load noisyt.mat
load noisys.mat

%
disp('times and signal values read in OK. No. of sample points is')
nt=length(t)

%make a plot of some of these
np=min(256,nt);nt2=nt/2;np2=np/2;
disp('Press RETURN for a plot of first so many values')
pause
xp=t(1:np); yp=s(1:np);
plot(xp, yp,'g+')
pause

% work out time scale and hence frequency scale
dt=t(2)-t(1); % the time spacing
f=(1:np2-1)/(dt*(np-1)); % the corrresponding frequency values
T=1 ./f; % a corresponding period scale for plotting

% Take the Fourier transform
yf=fft(yp,np); % take the Fourier transform
%
spect=yf(2:np2).*conj(yf(2:np2)); % this is the mod^2 of the
clg % frequency strength
disp(' Now plot against frequency....press RETURN')
pause

% put on a log/linear plot
semilogy(f,spect,'g+')
title(' Frequency spectrum')
xlabel('f /secs')
ylabel('P')
pause

% a log-log plot this time showing the periodicity
disp(' Now plot against period....press RETURN')
pause
loglog(T,spect,'r+')
title(' Periodicity spectrum')
xlabel('T secs')
ylabel('P')

fodesol obtains an approximate solution to the ODE

function fodesol(F,tmin,tmax,xmin,xmax)
%
%
% function fodesol(F,tmin,tmax,xmin,xmax)
%
% fodesol obtains an approximate solution to the ODE
% dx/dt = f(x,t) by the Runge-Kutta-Fehlberg
% method ODE23 supplied with MATLAB (2nd/3rd
% order version) and compares the solutions
% with a "grain" (or field) plot.
%
% F is a string containing the name of a real
% function f(x,t), supplied in an M-file.
% e.g. fnxt.m
% tmin,tmax and
% xmin,xmax are limits for the plot
%
% usage e.g.
% fodesol('fnxt',0, 2.5, -3 , 3)
%
% fodesol.m uses: odegr.m and mode23.m
%

% Initialization
disp('Here is a grain or field plot of dx/dt = f(x,t)')
disp('Press Enter to see it')
pause
more = 'y';
clg; hold off

% Grain plot
%axis([tmin tmax xmin xmax])
odegr(F,tmin,tmax,xmin,xmax)
hold on
axis(axis)
% Superimpose numerical solution through given point
while more == 'y'
t0 = input('Pick initial t-value: ');
x0 = input('Pick initial x-value: ');
disp('please wait...')
[xp,yp] = mode23(F,t0,tmax,x0);
% use a spline to make it look smoother
n=length(xp); x1=xp(1); x2=xp(n);
dx=(x2-x1)/50; xi=x1:dx:x2;
yi=spline(xp',yp',xi);
plot(xi,yi,'r'); pause
more = input('Do you want another solution? (y/n): ','s');
end
disp(' Initial and final points (t,x) of solution were')
l=length(yp);
tf=xp(l); xf=yp(l);
disp(' Initial Final')
disp([t0 tf])
disp([x0,xf])
hold off

Flu epidemic data

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%fludat.m
%
% Flu epidemic data
%
% Usage:
% fludat
%
% Loads in matrix alldat(ndata,2)
% and vectors xi (tdays), yi (I)
% which are its columns
% tdays is the time in days
% I is the number of infectives (ill boys) at time t
% The number of data points (ndata) is also returned
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

echo off;
hold off;

%Data goes in here
% tdays I
alldat=[
3 24
4 79
5 230
6 301
7 256
8 245
9 201
10 138
11 78
12 25
13 12
14 9];

% end of the data
% define a few convenient variables
tdays = alldat(:,1); xi=tdays;
I = alldat(:,2); yi=I;
disp('No. of data points:')
ndata = length(tdays);
disp(ndata)

%plot the data
disp(' Press RETURN for plot')
pause
%define limits of the plot (bottom LH, top RH)
xp=[0 16]; yp=[0 800];
plot(xp,yp,'.w')
hold on

plot(xi,yi,'*g')
title('Flu infections')
xlabel('days')
ylabel('I')
hold off

Estimates the convergence rate of a given sequence

function Rate_k = name(Root, Id);

% f_rate.m

%

% Estimates the convergence rate of a given sequence

% of scalars (Root=Vector) or vectors (Root=Matrix).

%

% Usage: K = f_rate ( Root, id)

% ==or== K = f_rate ( Root )

%

% where Root is of size (m x n), with n = dimension, and

% id=1 is to show all intermediate estimates.

% For usage 2, id=0;

% in all cases the very last estimate is shown.

% For example, to check the convergence of the iterates,

% (1, 2, 3, 4, 4.5, 5, 5.11)

% USE Root = [1 2 3 4 4.5 5 5.11]';

% kc_1 = f_rate(Root,1) %% Giving 1.3914 ... 2.8614

% kc_2 = f_rate(Root) %% Giving 1.7101 and 2.8614



if (nargin < 1), help f_rate; return, end

if nargin==1, Id=0; end

[m n] = size(Root) ; % Check if m is more than 4

if (m <>= 4), Root=Root'; [m n]=size(Root); end

abort_id = 0;

x_exact = Root(m,:) ; % Take as exact solution "x"

keep = [] ; Rate_k = []; i=3;

if m >= 4

while (i <> 3 or m = 3

d_n_m = norm(Root(i-2,:)-x_exact) ; % d_(n-1)

d_n = norm(Root(i-1,:)-x_exact) ; % d_n

d_n_p = norm(Root(i,:) -x_exact) ; % d_(n+1)

if (d_n_m < eps) d_n_m = eps*7; end % In case of 0

if (d_n < eps) d_n = eps*7; end % In case of 0

if (d_n_m < eps) d_n_m = eps*7; end % In case of 0

if (d_n < eps) d_n = eps*7; end % In case of 0

if (d_n_p < eps)

d_n_p = eps*7 ; % In case of 0

i = m - 1; abort_id = 1; % No need to continue

end

top = d_n_p / d_n ; % numerator

bottom = d_n / d_n_m ; % denominator

if (top == 1) top = 1 + eps*7 ; end % Safeguard for Log

if (bottom == 1) bottom = 1 + eps*7 ; end

rate_k = log(top) / log(bottom) ; % Estimate of the order

if (abort_id == 0 & rate_k > 0), keep = [keep; rate_k]; end

i = i + 1;

end % End WHILE-loop

disp(' ');

if Id == 1

Rate_k = keep ;

disp(' (Individual estimates for rate)') ;

else

m = length(keep); m=keep(m);

fprintf('Average Rate =%5.1f\n',mean(keep));

fprintf('Last estimate=%5.1f\n',m);

Rate_k = [mean(keep) m];

end % End for IF

else

fprintf('Rate_k is not estimated, because size = %d\n',m)

disp('Not enough iterates have been supplied for f_rate.m');

help f_rate ;

end; % END the main IF m>4

M-file to plot the Fibonacci numbers f(i) less than 1000

%M-file to plot the Fibonacci numbers f(i) less than 1000.
%The Fibonacci numbers are defined by f(1)=f(2)=1
%and f(i+2)=f(i+1)+f(i) for i=1,2,3,...
%They are calculated and then plotted with i along the
%horizontal axis and f(i) along the vertical axis.

f=[1 1];
k=1;
while f(k) < 1000
f(k+2)=f(k+1)+f(k);
k=k+1;
end
f
plot(f)