Не запускается пример из курса

clear all, clc, close all
%%%% Расчет площади эллипса методом Монте-Карло
%%% Параметры эллипса
a = 6;
b = 4;
x0 = -2;
y0 = 3;
 
%%% Построение эллипса
t = 0 : min([a, b])/500 : 2*pi;
x = a*cos(t) + x0;
y = b*sin(t) + y0;
line(x,y, 'linew', 2, 'color', [1,0,0])
%%% Центр пересечения осей эллипса
line(x0,y0,'marker','o','markerfacecolor','r','color','r' )
 
%%% Построение прямоугольника, в который вписан эллипс
line([x0-a, x0+a],[y0+b,y0+b], 'lines','-.')
line([x0-a, x0+a],[y0-b,y0-b], 'lines','-.')
line([x0-a, x0-a],[y0-b,y0+b], 'lines','-.')
line([x0+a, x0+a],[y0-b,y0+b], 'lines','-.')
 
%%% Гененрирование случайных чисел для области D* - области прямоугольника
N = 500; %%% число испытаний
%%% Генерация чисел по горизонтальной стороне прямоугольника
rx = (x0-a) + (x0+a - (x0-a))*rand(N,1);
%%% Генерация чисел по вертикальной стороне прямоугольника
ry = (y0-b) + (y0+b - (y0-b))*rand(N,1);
 
%%%% Заполнение случайными числами прямоугольника с эллипсом
for J = 1 : N
line(rx(J),ry(J),'marker','o','markersize',2,'color', 'k',...
 'markerfacecolor', 'k' )
end
 
%%% Подсчет количества случайных чисел, попавших в эллипс
m = 0;
for J = 1 : N
    if ((rx(J) - x0)/a)^2 + ((ry(J) - y0)/b)^2 <= 1
        m = m + 1;
    end
end
 
%%%%% Расчет площади прямоугольника
S = ((x0+a) - (x0-a))*((y0+b)-(y0-b));
 
%%% Расчет площади эллипса методом Монте-Карло
S0 = m/N*S
 
%%% Проверка расчета 
Scontrol = pi*a*b
 
%%% Заголовок для диаграммы
str = sprintf('%s Оси эллипса a = %g, b = %g, координаты центра x_0 = %g,  y_0 = %g.%s%g%s%g. %s%g', '\bf',a, b,x0,y0, ... 
'\newline Теоретическая площадь эллипса S = ', Scontrol, ...
'\newline Площадь эллипса по методу Монте-Карло S0 = ', S0, ...
 'Число испытаний N = ', N);
title(str)
grid off %% выключение сетки на диаграмме
xlabel('\bf - - - - - - - X  - - - - - - - ')
ylabel('\bf - - - - - - - Y  - - - - - - - ')
%%% Установление размеров и свойств графического окна
gfig = get(0, 'screensize');
set(gcf, 'color', 'w', 'position', [gfig(1) + 100, gfig(2) + 100, gfig(3)*0.8, gfig(4)*0.7])
%%% Ограничения по осям
xlim([x0-a-a/10, x0+a+a/10])
ylim([y0-b-b/10, y0+b+b/10])
axis equal