MATLAB stands for MATrix LABoratory. It is a high-level numerical computing environment and programming language by MathWorks. The .m extension refers to plain text files containing MATLAB code — either scripts or function definitions, collectively called M-files or M-scripts.

MATLAB excels at matrix operations, signal processing, control systems, simulation, and rapid algorithm prototyping.

% Script: myScript.m — variables go to base workspace
x = 10;
y = x ^ 2;

% Function: mySquare.m
function result = mySquare(x)
    result = x ^ 2;
end

Append a semicolon ; to suppress output. Without it, MATLAB echoes the result to the Command Window.

a = 5       % prints: a = 5
b = 10;     % suppressed
c = a + b;  % suppressed

• double — default 64-bit float (IEEE 754)
• single — 32-bit float
• int8/16/32/64, uint8/16/32/64 — signed/unsigned integers
• logical — boolean (true/false)
• char — character array (legacy strings)
• string — string array (R2016b+, recommended)
• cell — heterogeneous container
• struct — named-field data container

x = 3.14;           % double (default)
n = int32(100);     % 32-bit integer
flag = true;         % logical
s = "Hello";         % string
class(x)             % 'double'

The colon operator creates ranges and is one of MATLAB’s most powerful features — used for ranges, array indexing, and slicing.

v = 1:5           % [1 2 3 4 5]
v = 0:0.5:2       % [0 0.5 1.0 1.5 2.0]  (start:step:end)
v = 10:-2:2      % [10 8 6 4 2]

A = magic(4);
A(:, 2)            % all rows, column 2
A(1, :)            % row 1, all columns
A(:)               % all elements as column vector

• = — assignment operator
• == — equality comparison (returns logical array)
• ~= — not-equal comparison

x = 5;             % assignment
x == 5             % true  (1)
x ~= 3             % true  (1)
[1 2 3] == 2      % [0 1 0] — element-wise

clear              % clear all workspace variables
clear x y          % clear specific variables
clc                % clear Command Window display
close all         % close all figure windows
clearvars -except x  % clear all except x

ans is MATLAB’s automatic variable that stores the result of the most recent unassigned expression.

3 + 4           % ans = 7
ans * 2         % ans = 14

* is matrix multiplication; .* is element-wise multiplication. Same distinction applies to / vs ./ and ^ vs .^.

A = [1 2; 3 4]; B = [5 6; 7 8];
A * B     % matrix product: [19 22; 43 50]
A .* B    % element-wise:   [5 12; 21 32]

v = [1 2 3];
v .^ 2    % [1 4 9]
v ^ 2     % ERROR — v is not square

disp('Hello World')               % simple display
fprintf('Value: %.2f\n', 3.14159) % formatted: "Value: 3.14"
sprintf('x = %d', 42)             % returns formatted string

• eps — machine epsilon (~2.22e-16); smallest representable difference from 1.0
• inf — positive infinity (result of 1/0)
• nan — Not a Number (result of 0/0 or inf-inf)
• pi — π ≈ 3.14159265358979

1/0          % Inf
0/0          % NaN
isnan(0/0)  % true
isinf(1/0)  % true

help fft           % shows help text in Command Window
doc fft            % opens full documentation browser
lookfor 'inverse'  % searches all help text for keyword

who lists variable names; whos shows detailed info: size, bytes, and class.

x = rand(100);
whos x
%  Name    Size          Bytes  Class
%  x       100x100       80000  double

% Single-line comment

%{
  Block comment
  (multiple lines)
%}

x = 5; % inline comment after code

result = 1 + 2 + 3 + ...
         4 + 5 + 6;   % 21

A = [1 2 3; ...
     4 5 6];

Three dots ... (ellipsis) are MATLAB’s line continuation operator.

A = [1 2 3; 4 5 6; 7 8 9]   % 3×3 (space/comma = col, ; = row)
B = zeros(3, 4)                   % 3×4 all-zeros
C = ones(2)                       % 2×2 all-ones
D = eye(4)                        % 4×4 identity
E = rand(3)                       % 3×3 uniform random [0,1]
F = randn(3)                      % 3×3 standard normal

r = [1 2 3]       % row vector    (1×3)
c = [1;2;3]       % column vector (3×1)

r'                 % conjugate transpose (use for complex data)
r.'                % non-conjugate transpose (safe for real data)
transpose(r)      % equivalent to r.'

A = [10 20 30; 40 50 60; 70 80 90];
A(2, 3)           % element row 2, col 3 → 60
A(1, :)           % entire row 1    → [10 20 30]
A(:, 2)           % entire col 2    → [20;50;80]
A(1:2, 2:3)     % submatrix rows 1-2, cols 2-3
A(end, :)         % last row
A(end-1:end, :) % last two rows

MATLAB uses 1-based indexing. The first element is at index 1. This is a common source of off-by-one errors for programmers from C or Python.

v = [10 20 30];
v(1)    % 10 — first element
v(end)  % 30 — last element
v(0)    % ERROR: index must be a positive integer

A = rand(3, 5);
size(A)          % [3 5]
size(A, 1)       % 3  (rows)
size(A, 2)       % 5  (columns)
length(A)        % 5  (largest dimension)
numel(A)         % 15 (total elements)
ndims(A)         % 2  (number of dimensions)

linspace(0, 1, 5)   % [0 0.25 0.5 0.75 1.0] — exactly 5 points
0:0.25:1             % same, but relies on step arithmetic

linspace(a,b,n) guarantees exactly n points — safer than the colon operator when floating-point rounding could matter.

A = [1 2 3; 4 5 6];    % 2×3
reshape(A, 3, 2)          % 3×2 (column-major)
reshape(A, 1, [])         % 1×6 row vector (auto-infer)
A(:)                       % 6×1 column vector

A = [1 2]; B = [3 4];
[A B]           % horizontal: [1 2 3 4]
[A; B]          % vertical:   [1 2; 3 4]
cat(3, A, B)   % along 3rd dimension

v = [3 -1 5 -2 4];
mask = v > 0;         % [1 0 1 0 1] logical array
v(mask)               % [3 5 4] — select positives
v(v > 0)             % same, inline
v(v < 0) = 0;        % set negatives to zero

A = magic(4);
A(A > 10)            % all elements > 10 as column vector

v = [5 0 3 0 7];
find(v)               % [1 3 5] — indices of nonzero elements
find(v > 3)          % [1 5]   — indices where v > 3
find(v, 2)           % [1 3]   — first 2 nonzero indices
[r, c] = find(A > 5) % row and column indices for matrices

A = [2 1; 5 3];
inv(A)            % matrix inverse
det(A)            % determinant → 1
rank(A)           % rank → 2
trace(A)          % sum of diagonal → 5
[V, D] = eig(A)  % V = eigenvectors, D = eigenvalue diagonal
[U, S, V] = svd(A) % singular value decomposition

v = [3 1 4 1 5 9];
sort(v)                    % [1 1 3 4 5 9] ascending
sort(v, 'descend')        % [9 5 4 3 1 1]
[sorted, idx] = sort(v)   % idx = original positions
sortrows(A, 2)            % sort matrix rows by column 2

Sparse matrices store only non-zero elements, saving memory and computation when most entries are zero (e.g., adjacency matrices, FEM stiffness matrices).

S = sparse([1 2 3], [1 2 3], [10 20 30], 5, 5)
% 5×5 sparse diagonal matrix
full(S)           % convert to full matrix
nnz(S)            % number of non-zero elements
spy(S)            % visualise sparsity pattern

x = 75;
if x >= 90
    grade = 'A';
elseif x >= 80
    grade = 'B';
elseif x >= 70
    grade = 'C';
else
    grade = 'F';
end

Note: MATLAB uses elseif (one word), not else if.

for i = 1:5
    fprintf('i = %d\n', i);
end

% Loop over columns of a matrix:
A = [1 2 3; 4 5 6];
for col = A       % col is each column: [1;4], [2;5], [3;6]
    disp(col);
end

x = 1;
while x < 100
    x = x * 2;
end
disp(x)    % 128

for i = 1:10
    if i == 3
        continue    % skip this iteration
    end
    if i == 7
        break       % exit loop entirely
    end
    fprintf('%d ', i);
end
% Output: 1 2 4 5 6

day = 'Mon';
switch day
    case 'Mon'
        disp('Monday');
    case {'Sat', 'Sun'}   % multiple values in one case
        disp('Weekend');
    otherwise
        disp('Weekday');
end

Unlike C, MATLAB’s switch does not fall through — no break needed.

Without pre-allocation, MATLAB dynamically grows the array each iteration by copying the entire array — O(n²) time. Pre-allocation is O(n).

% BAD — grows every iteration (slow)
result = [];
for i = 1:10000
    result(i) = i^2;
end

% GOOD — pre-allocated (fast)
result = zeros(1, 10000);
for i = 1:10000
    result(i) = i^2;
end

% Loop version (slow)
x = 1:1000; y = zeros(1,1000);
for i = 1:1000
    y(i) = sin(x(i)) * exp(-x(i)/10);
end

% Vectorized (10-100x faster)
x = 1:1000;
y = sin(x) .* exp(-x/10);

Vectorisation leverages MATLAB’s optimised BLAS/LAPACK routines and JIT compilation.

x = [1 2 3 4];
y = arrayfun(@(n) n^2 + 1, x)   % [2 5 10 17]

Use when the operation can’t be vectorised (e.g. function only accepts scalars). For performance-critical code, a proper vectorised approach or pre-allocated loop is usually faster.

words = {'hello', 'world', 'matlab'};
lens  = cellfun(@length, words)      % [5 5 6]
upper = cellfun(@upper, words, 'UniformOutput', false)
% {'HELLO','WORLD','MATLAB'} — UniformOutput=false for cell results

try
    x = inv([1 2; 1 2]);   % singular matrix
    error('Custom error: %s', 'something failed');
catch e
    fprintf('Caught: %s\n', e.message);
    fprintf('ID:     %s\n', e.identifier);
end

function [mn, mx, rng] = stats(v)
    mn  = min(v);
    mx  = max(v);
    rng = mx - mn;
end

[a, b, c] = stats([3 1 4 1 5]);  % a=1, b=5, c=4
[a, ~, c] = stats([3 1 5]);       % ~ discards 2nd output

function result = myFunc(x, y, z)
    if nargin < 3, z = 0; end
    if nargin < 2, y = 1; end
    result = x + y + z;
end

myFunc(5)          % 6  (x=5, y=1, z=0)
myFunc(5, 2)       % 7
myFunc(5, 2, 3)    % 10

square = @(x) x.^2;
square(5)               % 25
square([1 2 3])         % [1 4 9]

add = @(x, y) x + y;
add(3, 4)               % 7

% Variables are captured at creation time:
a = 10;
f = @(x) x + a;
a = 99;
f(5)                     % 15 (captured a=10, not 99)

fh = @sin;              % handle to built-in
fh(pi/2)               % 1.0
fplot(@cos, [0 2*pi])

% Pass as argument to numerical solvers:
root = fzero(@(x) x^2 - 2, 1)   % √2 ≈ 1.4142

Local functions are helpers defined after the main function in the same .m file — private to that file. Nested functions are defined inside another function and share its workspace.

function y = main(x)
    y = helper(x) * 2;
end

function r = helper(x)   % local function
    r = x + 1;
end

function result = sumAll(varargin)
    result = 0;
    for k = 1:nargin
        result = result + varargin{k};
    end
end

sumAll(1, 2, 3, 4)   % 10

function f = factorial_r(n)
    if n <= 1
        f = 1;
    else
        f = n * factorial_r(n - 1);
    end
end

factorial_r(5)   % 120

function result = myFunc(x, varargin)
    p = inputParser();
    addRequired(p, 'x', @isnumeric);
    addOptional(p, 'scale', 1, @isnumeric);
    addParameter(p, 'method', 'linear', @ischar);
    parse(p, x, varargin{:});
    result = p.Results.x * p.Results.scale;
end

myFunc(5, 2, 'method', 'cubic')  % result = 10

function y = safeSqrt(x)
    if x < 0
        y = NaN;
        return;          % exit function early
    end
    y = sqrt(x);
end

function count = callCounter()
    persistent n;
    if isempty(n), n = 0; end
    n = n + 1;
    count = n;
end
callCounter()  % 1
callCounter()  % 2

% fzero expects a function of one variable
% but our equation has extra parameter 'a':
a = 3;
f = @(x) x^2 - a;      % 'a' captured in closure
root = fzero(f, 1)      % √3 ≈ 1.7321

f = @(x) x.^2;
g = @(x) x + 1;
h = @(x) f(g(x));    % h(x) = (x+1)^2
h(3)                 % 16

profile on
myScript();
profile viewer    % opens interactive HTML profiler
profile off

The profiler shows time per function and per line — essential before any optimisation effort.

function y = myFunc(x, n)
    validateattributes(x, {'numeric'}, {'vector','nonempty'}, 'myFunc', 'x');
    validateattributes(n, {'numeric'}, {'scalar','positive','integer'});
    y = x(n);
end

Generates clear error messages like: “Expected input x to be a vector, but its size is [3 3].”

c = 'hello'; c(1)  % 'h' — character
s = "hello"; s(1)  % "hello" — whole string
"foo" + " bar"     % "foo bar"
"foo" == "foo"     % true

A cell array holds elements of different types or sizes. Use it for heterogeneous data.

C = {42, 'hello', [1 2 3], true};
C{1}       % 42      — curly braces extract content
C(1)       % {42}    — round brackets return sub-cell
C{3}(2)   % 2       — index into extracted content

s = "Hello World";
strrep(s, "World", "MATLAB")   % "Hello MATLAB"
contains(s, "World")            % true
startsWith(s, "Hello")          % true
endsWith(s, "World")            % true
strfind(s, "o")                 % [5 8] — positions (1-based)

s = "a,b,c,d";
parts  = strsplit(s, ",")        % {"a","b","c","d"}
joined = strjoin(parts, " - ")  % "a - b - c - d"

num2str(3.14)        % '3.14' (char)
str2double("3.14")   % 3.14 — recommended, works on string arrays
str2num('42')        % 42   — uses eval internally (avoid)
int2str(7.9)         % '8'  — rounds to nearest integer

s = 'Order #1234 placed on 2024-06-15';

regexp(s, '\d+', 'match')     % {'1234','2024','06','15'}

% Named tokens:
tok = regexp(s, '(?<year>\d{4})-(?<month>\d{2})', 'names');
tok.year     % '2024'

% Replace:
regexprep(s, '\d+', 'N')      % 'Order #N placed on N-N-N'

person.name   = 'Alice';
person.age    = 30;
person.scores = [85 92 78];

% Or with struct():
person = struct('name','Alice', 'age',30);

fieldnames(person)          % {'name','age','scores'}
rmfield(person, 'age')     % remove field
isfield(person, 'name')   % true

people(1) = struct('name','Alice','age',30);
people(2) = struct('name','Bob',  'age',25);

people(1).name      % 'Alice'
[people.age]        % [30 25] — extract field from all elements

m = containers.Map({'apple','banana'}, {1.5, 0.75});
m('apple')           % 1.5
m('cherry') = 2.0;  % add key
isKey(m, 'banana')   % true
keys(m)              % {'apple','banana','cherry'}
values(m)            % {1.5, 0.75, 2.0}

% For char arrays — use strcmp (not ==):
strcmp('hello', 'hello')          % true
strcmpi('Hello', 'hello')         % true (case-insensitive)

% For string arrays — == works:
"hello" == "hello"               % true
isequal("foo", "foo")            % true (all types)

% Read
T = readtable('data.csv');          % as table (recommended)
M = readmatrix('data.csv');        % as numeric matrix

% Write
writetable(T, 'output.csv');
writematrix(M, 'output.csv');

x = 42; A = rand(3);
save('mydata.mat')              % save all variables
save('mydata.mat', 'x', 'A')  % save specific variables
save('big.mat', '-v7.3')       % HDF5 format for >2 GB files

load('mydata.mat')              % loads into workspace
s = load('mydata.mat')          % loads into struct: s.x, s.A

fid = fopen('myfile.txt', 'r');
if fid == -1, error('Cannot open file'); end

lines = {};
while ~feof(fid)
    line = fgetl(fid);
    if ischar(line)
        lines{end+1} = line;
    end
end
fclose(fid);

fgetl returns a char array for each line, or -1 at end-of-file. Always call fclose to release the file handle.

T = readtable('data.xlsx');
T = readtable('data.xlsx', 'Sheet', 'Sales', 'Range', 'A1:E100');
writetable(T, 'output.xlsx', 'Sheet', 'Results');

fid = fopen('results.txt', 'w');
fprintf(fid, '%-10s %8s %8s\n', 'Name', 'Score', 'Grade');
fprintf(fid, '%-10s %8.2f %8s\n', 'Alice', 92.5, 'A');
fclose(fid);

files = dir('data/*.csv');
for k = 1:numel(files)
    fname = fullfile(files(k).folder, files(k).name);
    T = readtable(fname);
    fprintf('Loaded %s: %d rows\n', files(k).name, height(T));
end

T = readtable('sales.csv');

T.Revenue                              % access a column as vector
T(T.Revenue > 1000, :)              % filter rows
T(strcmp(T.Region, 'APAC'), :)      % filter by string column
T.Profit = T.Revenue - T.Cost;        % add computed column
G = groupsummary(T, 'Region', 'sum', 'Revenue'); % group-by

% Read JSON from file (R2016b+)
data = jsondecode(fileread('data.json'));

% Write JSON to file
s = struct('name', 'Alice', 'scores', [90 85]);
json = jsonencode(s);
fid = fopen('out.json', 'w');
fprintf(fid, '%s', json);
fclose(fid);

x = linspace(0, 2*pi, 200);
figure;
plot(x, sin(x), 'b-', 'LineWidth', 2);
hold on
plot(x, cos(x), 'r--', 'LineWidth', 2);
hold off
xlabel('x (rad)'); ylabel('Amplitude');
title('Sine and Cosine');
legend('sin(x)', 'cos(x)');
grid on;

x = 0:0.1:2*pi;
figure('Position', [100 100 800 400]);

subplot(1, 2, 1);    % 1 row, 2 cols, panel 1
plot(x, sin(x)); title('sin');

subplot(1, 2, 2);    % panel 2
plot(x, cos(x)); title('cos');

plot(x, y)          % line plot
scatter(x, y)       % scatter plot
bar(x, y)           % bar chart
histogram(data)     % histogram
boxplot(data)       % box-and-whisker
pie(values)         % pie chart
surf(X, Y, Z)       % 3D surface
contour(X, Y, Z)   % contour lines
imagesc(A)          % colour-mapped matrix
heatmap(T, r, c)   % heatmap from table

fig = gcf;    % get current figure handle
saveas(fig, 'myplot.png')
saveas(fig, 'myplot.pdf')
exportgraphics(fig, 'myplot.pdf', 'ContentType', 'vector')
print(fig, 'myplot', '-dpng', '-r300')  % 300 DPI PNG

[X, Y] = meshgrid(-3:0.1:3);
Z = sin(sqrt(X.^2 + Y.^2));

figure;
surf(X, Y, Z);
shading interp;      % smooth shading (no grid lines)
colormap parula;
colorbar;
xlabel('X'); ylabel('Y'); zlabel('Z');

By default each plot call clears the axes. hold on preserves existing plots so subsequent calls overlay on the same axes. hold off restores the default replace behaviour.

plot(x, sin(x));
hold on
plot(x, cos(x));   % overlaid on same axes
hold off

xlim([0 10]);         % x-axis range
ylim([-1 1]);        % y-axis range
axis([0 10 -1 1]);  % [xmin xmax ymin ymax]
axis equal;           % equal aspect ratio
axis tight;           % fit tightly to data
axis off;             % hide axes entirely

x = 0:0.1:2*pi;
plot(x, sin(x));

text(pi/2, 1.05, '\leftarrow Peak', 'FontSize', 12)
annotation('arrow', [0.4 0.5], [0.8 0.9])  % normalised coords
legend('sin(x)', 'Location', 'northeast');

A = [2 1; 5 3];
b = [8; 19];

x = A \ b           % backslash — recommended [1; 6]
x = inv(A) * b     % avoid (slower, less numerically stable)
x = linsolve(A, b) % with structural hints (opts)

% x^3 - 6x^2 + 11x - 6 = 0  (roots at 1, 2, 3)
p = [1 -6 11 -6];      % coefficients, highest power first
roots(p)               % [3; 2; 1]

polyval(p, 2)         % evaluate polynomial at x=2 → 0
polyfit(x, y, 3)      % fit degree-3 polynomial to data

f = @(x) x^3 - x - 2;
root = fzero(f, 1)           % initial guess → 1.5214
root = fzero(f, [1 2])       % bracket where sign changes

opts = optimset('TolX', 1e-12, 'Display', 'off');
fzero(f, 1, opts)

% Single variable — bounded interval:
f = @(x) (x-3).^2 + 1;
[xmin, fval] = fminbnd(f, 0, 10)   % xmin=3, fval=1

% Multi-variable (derivative-free):
g = @(v) (v(1)-1)^2 + (v(2)-2)^2;
[xmin, fval] = fminsearch(g, [0, 0])

% With Optimization Toolbox (gradient-based):
[xmin, fval] = fminunc(g, [0 0])

f = @(x) sin(x) ./ x;
I = integral(f, 0.01, pi)      % adaptive quadrature

% 2D integral:
g = @(x,y) x .* y.^2;
I2 = integral2(g, 0, 1, 0, 1)

% Trapezoidal rule on sampled data:
x = 0:0.01:1;
trapz(x, sin(x))

% dy/dt = -2y, y(0) = 1 → exact solution: e^(-2t)
f = @(t, y) -2 * y;
[t, y] = ode45(f, [0 5], 1);

plot(t, y, t, exp(-2*t), 'r--');
legend('ode45', 'exact');

Use ode23 for lower accuracy, ode15s for stiff problems, ode113 for high accuracy smooth solutions.

Fs = 1000; N = 1024;
t  = (0:N-1)/Fs;
x  = sin(2*pi*50*t) + sin(2*pi*120*t);  % 50 + 120 Hz tones

X = fft(x);
f = (0:N-1) * (Fs/N);
plot(f(1:N/2), abs(X(1:N/2)));    % one-sided spectrum
xlabel('Frequency (Hz)');

x = [1 2 3 4 5];
y = [2.1 4.0 5.9 8.1 10.2];
p = polyfit(x, y, 1)        % linear fit: p=[slope intercept]
y_fit = polyval(p, x)       % evaluate at x

% Nonlinear (Curve Fitting Toolbox):
f = fittype('a*exp(-b*x)');
[cf, gof] = fit(x', y', f, 'StartPoint', [1 0.5]);

x = [0 1 2 3 4];
y = [0 1 4 9 16];

interp1(x, y, 1.5)                   % 2.5  (linear, default)
interp1(x, y, 1.5, 'spline')         % 2.25 (cubic spline)
interp1(x, y, 1.5, 'pchip')          % shape-preserving cubic
interp1(x, y, 1.5, 'nearest')        % nearest neighbour
interp2(X, Y, Z, xi, yi)             % 2D grid interpolation

data = randn(1, 1000);
mean(data)              % sample mean
std(data)               % standard deviation (n-1)
var(data)               % variance
median(data)            % median
prctile(data, [25 75]) % 25th and 75th percentile
corrcoef(x, y)          % correlation matrix
cov(x, y)              % covariance matrix

% File: BankAccount.m
classdef BankAccount
    properties
        Owner  string
        Balance double = 0
    end
    methods
        function obj = BankAccount(owner, bal)
            obj.Owner   = owner;
            obj.Balance = bal;
        end
        function obj = deposit(obj, amount)
            obj.Balance = obj.Balance + amount;
        end
    end
end

acc = BankAccount("Alice", 1000);
acc = acc.deposit(500);    % acc.Balance = 1500

classdef Counter < handle
    properties, Value = 0; end
    methods
        function increment(obj)   % no return value needed
            obj.Value = obj.Value + 1;
        end
    end
end

c1 = Counter; c2 = c1;   % c2 points to SAME object
c1.increment();
c2.Value   % 1 — c2 sees the change

classdef Animal
    properties, Name string; end
    methods
        function speak(obj)
            disp([obj.Name ' makes a sound']);
        end
    end
end

classdef Dog < Animal      % inherits from Animal
    methods
        function speak(obj)  % override method
            disp([obj.Name ' says: Woof!']);
        end
    end
end

d = Dog; d.Name = "Rex"; d.speak()  % Rex says: Woof!

classdef Sensor < handle
    events
        DataReady
    end
    methods
        function read(obj)
            notify(obj, 'DataReady');   % fire the event
        end
    end
end

s = Sensor;
addlistener(s, 'DataReady', @(~,~) disp('Data received!'));
s.read();   % triggers listener → prints 'Data received!'

classdef Vec2
    properties, x, y; end
    methods
        function obj = Vec2(x,y), obj.x=x; obj.y=y; end
        function r = plus(a, b)     % overloads the + operator
            r = Vec2(a.x+b.x, a.y+b.y);
        end
        function disp(obj)
            fprintf('(%g, %g)\n', obj.x, obj.y);
        end
    end
end

Vec2(1,2) + Vec2(3,4)   % (4, 6)

MATLAB uses copy-on-write: assignment is cheap until a modification triggers a deep copy. Key concerns:

• Growing arrays in loops without pre-allocation (O(n²) copies)
• Passing large arrays to functions that modify them (triggers copy)
• Accumulating large results without clearing intermediate variables

whos           % inspect variable sizes and bytes
clear largeVar  % release memory explicitly
memory         % total MATLAB memory usage (Windows)

syms x
f = x^3 - 6*x^2 + 11*x - 6;

factor(f)         % (x-1)*(x-2)*(x-3)
diff(f, x)        % 3x^2 - 12x + 11
int(f, x)         % indefinite integral
int(f, x, 0, 1)  % definite integral from 0 to 1
solve(f == 0, x)  % [1; 2; 3]
limit(f/x, x, 0) % symbolic limit as x→0

% parfor — parallel for (Parallel Computing Toolbox)
parpool(4);     % open pool of 4 workers

result = zeros(1, 1000);
parfor i = 1:1000
    result(i) = someHeavyFn(i);  % iterations run in parallel
end
delete(gcp);    % close pool

% GPU acceleration:
A_gpu = gpuArray(rand(1000));
B_gpu = A_gpu * A_gpu;           % computed on GPU
B     = gather(B_gpu);          % transfer result back to CPU

MEX files are compiled C/C++/Fortran routines that behave like MATLAB functions — used when vectorisation can’t achieve required speed (e.g., tight loops with data dependencies).

% 1. Write myMex.c with mexFunction() as the entry point
% 2. Compile in MATLAB:
mex myMex.c        % creates myMex.mex64 (Windows)

% 3. Call like any built-in:
result = myMex(inputData);

Simulink is a graphical block-diagram environment for modelling and simulating dynamic systems (control, signal processing, physical systems). M-scripts interact with Simulink to:

• Set model parameters: set_param('model/Block','Gain','5')
• Run simulations programmatically: simOut = sim('myModel')
• Post-process and plot logged signals from simOut
• Embed MATLAB code directly via MATLAB Function blocks

A = [1 2]; B = [1 2];
A == B            % [1 1] — element-wise, returns logical array
isequal(A, B)    % true  — whole-value equality, scalar result

% NaN handling:
NaN == NaN        % false (IEEE 754 standard)
isequal(NaN, NaN)  % false
isequaln(NaN, NaN) % true  (NaN-aware comparison)

Since R2016b, MATLAB supports implicit expansion: element-wise operations between arrays of compatible (non-identical) sizes expand automatically — like NumPy broadcasting.

A = [1;2;3];       % 3×1 column
B = [10 20 30];    % 1×3 row
A + B
% Result: 3×3 matrix
%  [11 21 31
%   12 22 32
%   13 23 33]

% Pre-R2016b equivalent using bsxfun:
bsxfun(@plus, A, B)   % still works in all versions

% Click the dash to the left of a line in the Editor, or:
dbstop in myFunction at 25       % line 25
dbstop in myFunction if x < 0   % conditional breakpoint
dbstop if error                  % break on any error
dbclear all                      % remove all breakpoints

% At the K>> debug prompt:
%  dbstep    — step one line
%  dbstep in — step into a function call
%  dbcont    — continue to next breakpoint
%  dbquit    — exit debug mode

% tic/toc — wall-clock elapsed time
tic;
myHeavyFunction();
t = toc;
fprintf('Elapsed: %.4f s\n', t);

% timeit — more accurate; averages many runs and accounts for warm-up
t = timeit(@() myHeavyFunction());

MATLAB’s Code Analyzer highlights potential bugs, style issues, and performance warnings in the Editor (orange/red underlines). Run programmatically with:

checkcode('myScript.m')
checkcode('myScript.m', '-fullpath', '-id')  % include message IDs

Common warnings: array growing in loop, unused variables, using i/j as loop variables, missing semicolons, deprecated functions.

In MATLAB, i and j represent the imaginary unit (√−1). Using them as loop counters shadows the built-in, causing silent errors in complex arithmetic code.

for i = 1:5   % BAD — shadows imaginary unit
for k = 1:5   % GOOD — use k, ii, idx, etc.

z = 3 + 4*1i   % always use numeric literal 1i for imaginary

% File: TestMySquare.m
classdef TestMySquare < matlab.unittest.TestCase
    methods(Test)
        function testPositive(tc)
            tc.assertEqual(mySquare(3), 9)
        end
        function testNegative(tc)
            tc.assertEqual(mySquare(-4), 16)
        end
        function testTolerance(tc)
            tc.assertAlmostEqual(mySquare(1.5), 2.25, 'AbsTol', 1e-10)
        end
    end
end

% Run: runtests('TestMySquare')

% Compile myFunction to C:
codegen myFunction -args {zeros(1,100,'double')} -report

% Key constraints for codegen-compatible code:
% - Types must be statically inferable (no dynamic typing)
% - No cell arrays (use fixed-size arrays)
% - No eval/feval
% - Struct fields must be fixed at compile time

• Signal Processing Toolbox — filters, spectral analysis, wavelets
• Control System Toolbox — Bode plots, PID design, transfer functions
• Statistics & Machine Learning Toolbox — regression, classification, clustering
• Deep Learning Toolbox — CNNs, LSTMs, training loops, GPU support
• Optimization Toolbox — linear/nonlinear/integer programming
• Image Processing Toolbox — filtering, morphology, segmentation
• Parallel Computing Toolbox — parfor, GPU arrays, cluster jobs
• Symbolic Math Toolbox — symbolic calculus, algebra, equation solving

groups = [1;2;1;3;2;1];
vals   = [10;20;30;40;50;60];

accumarray(groups, vals)          % sum per group: [100; 70; 40]
accumarray(groups, vals, [], @max) % max per group: [60; 50; 40]
accumarray(groups, vals, [], @mean)% mean per group

accumarray is highly optimised — far faster than loop-based group aggregation on large datasets.

str2double("3.14")                   % 3.14
str2double(["1.5","2.5","bad"])    % [1.5, 2.5, NaN]

Use App Designer (R2016a+, recommended) which generates a class-based .mlapp file with auto-generated layout code. The legacy GUIDE tool is also available.

% Open App Designer:
appdesigner

% Programmatic UI (no designer needed):
fig = uifigure('Name', 'My App', 'Position', [100 100 400 300]);
btn = uibutton(fig, 'push', ...
    'Text', 'Click Me', ...
    'Position', [150 130 100 40], ...
    'ButtonPushedFcn', @(~,~) msgbox('Hello!'));
lbl = uilabel(fig, 'Text', 'Status: Ready', 'Position', [120 80 200 30]);

LU decomposition factors a matrix A into a lower triangular L and upper triangular U (with optional permutation P), so that PA = LU. MATLAB uses it internally for backslash.

A = [2 1 1; 4 3 3; 8 7 9];
[L, U, P] = lu(A)     % PA = LU

% Solve Ax = b using LU (faster when solving many rhs):
b = [1; 2; 3];
y = L \ (P*b);
x = U \ y;    % equivalent to A\b but reuses factorisation

QR factorises A into an orthogonal Q and upper triangular R. Used for least-squares, eigenvalue algorithms, and solving rank-deficient systems.

[Q, R] = qr(A)          % full QR
[Q, R, E] = qr(A)       % with column pivoting

% Least-squares: min ||Ax - b||
x = R \ (Q' * b)         % equivalent to A\b for overdetermined systems

% Economy QR (only r columns of Q):
[Q, R] = qr(A, 0)

Cholesky factorises a symmetric positive-definite matrix as A = R’*R (R upper triangular). It is twice as fast as LU for such matrices.

A = [4 2; 2 3];       % symmetric positive-definite
R = chol(A)             % upper triangular R: A = R'*R
[R, p] = chol(A)       % p=0 means A is pos-def; p>0 means failure

% Solve Ax = b efficiently:
x = R \ (R' \ b);

The condition number measures how sensitive the solution of Ax = b is to perturbations. A large condition number means the system is ill-conditioned and small errors in data cause large errors in the solution.

A = hilb(8);           % notoriously ill-conditioned
cond(A)                 % 2-norm condition number (~1.5e10 for hilb(8))
condest(A)              % fast estimate for large sparse matrices
rcond(A)                % reciprocal condition estimate (1/cond)

% Rule of thumb: lose ~log10(cond) digits of accuracy

The backslash operator (\) is not a single algorithm — it inspects the matrix and picks the most efficient solver:

• Square, full-rank: LU with partial pivoting
• Symmetric positive-definite: Cholesky
• Triangular: forward/back substitution
• Overdetermined (m > n): QR least-squares
• Underdetermined (m < n): minimum-norm QR
• Sparse: UMFPACK or SuiteSparse solvers

A \ b     % MATLAB chooses optimally — always prefer over inv(A)*b

rand(3)           % uniform [0,1]
randn(3)          % standard normal N(0,1)
randi(10, 3)      % uniform integers 1..10

% Reproducible results — seed the RNG:
rng(42)            % set seed (Mersenne Twister by default)
rng('default')     % reset to MATLAB default

% From Statistics Toolbox distributions:
normrnd(5, 2, 100, 1)   % 100 samples N(5,4)
exprnd(3, 100, 1)        % exponential, mean 3
binornd(10, 0.5, 50, 1) % binomial B(10,0.5)

x = 0:0.01:2*pi;
y = sin(x);

% First-order finite differences:
dy = diff(y) ./ diff(x)   % forward difference (length n-1)

% gradient uses central differences (same length as y):
dy = gradient(y, x(2)-x(1))

% For functions, use symbolic diff or complex-step differentiation:
f = @(x) sin(x.^2);
h = 1e-8;
df = (f(x+h) - f(x-h)) / (2*h);  % central finite difference

% System: x^2 + y^2 = 4,  x*y = 1
F = @(v) [v(1)^2 + v(2)^2 - 4;
          v(1)*v(2) - 1];

x0 = [1; 1];     % initial guess
sol = fsolve(F, x0)   % requires Optimization Toolbox

% Check residual:
F(sol)   % should be close to [0; 0]

The Kronecker product tiles one matrix with scaled copies of another. Common in quantum computing, signal processing, and solving Sylvester equations.

A = [1 0; 0 1];
B = [1 2; 3 4];
K = kron(A, B)
% [1 2 0 0
%  3 4 0 0
%  0 0 1 2
%  0 0 3 4]

% 2D finite-difference Laplacian on n×n grid:
T = toeplitz([-2 1 zeros(1,n-2)]);
L = kron(eye(n), T) + kron(T, eye(n));

a = [1 2 3];
b = [0 1 0.5];

conv(a, b)             % linear convolution (length m+n-1)
conv(a, b, 'same')    % same length as a
conv(a, b, 'valid')   % only fully overlapping region

xcorr(a, b)            % cross-correlation (Signal Processing Toolbox)
xcorr(a)               % autocorrelation

% Fast convolution via FFT:
N = length(a) + length(b) - 1;
c = ifft(fft(a,N) .* fft(b,N));

% y'' + y = 0,  y(0)=0, y(pi)=0  (BVP, not IVP)
% Rewrite as system: y1'=y2, y2'=-y1
odefun = @(x,y) [y(2); -y(1)];
bcfun  = @(ya,yb) [ya(1); yb(1)];  % y(0)=0, y(pi)=0

xmesh = linspace(0, pi, 20);
solinit = bvpinit(xmesh, [0; 1]);  % initial guess
sol = bvp4c(odefun, bcfun, solinit);

x = linspace(0, pi, 200);
y = deval(sol, x);
plot(x, y(1,:));

% Solve Laplace equation ∇²u = 0 on unit disk
model = createpde();
geometryFromEdges(model, @circleg);
applyBoundaryCondition(model, 'dirichlet', 'Edge', 1:4, 'u', 0);
specifyCoefficients(model, 'm',0,'d',0,'c',1,'a',0,'f',1);
generateMesh(model, 'Hmax', 0.1);
results = solvepde(model);
pdeplot(model, 'XYData', results.NodalSolution);

A timetable is like a table but with a datetime row index — ideal for time-series data.

t = datetime(2024,1,1) + hours(0:23)';
temp = 20 + 5*randn(24,1);
TT = timetable(t, temp, 'VariableNames', {'Temperature'});

TT(1:6, :)                           % first 6 hours
TT(TT.Temperature > 22, :)           % filter by value
TT_h = retime(TT, 'daily', 'mean')  % resample to daily mean

A = table([1;2;3], {'a';'b';'c'}, 'VariableNames', {'id','name'});
B = table([1;2;4], [100;200;400], 'VariableNames', {'id','score'});

innerjoin(A, B, 'Keys', 'id')   % rows present in both (ids 1,2)
outerjoin(A, B, 'Keys', 'id', 'MergeKeys', true)  % all rows
join(A, B)                        % left join on matching keys

T = table({'A';'B';'A';'B';'A'}, [10;20;30;40;50], ...
    'VariableNames', {'Group','Value'});

G = groupsummary(T, 'Group', {'sum','mean','max'}, 'Value')
%  Group  GroupCount  sum_Value  mean_Value  max_Value
%  'A'    3           90         30          50
%  'B'    2           60         30          40

This is equivalent to SQL: SELECT Group, COUNT(*), SUM(Value), AVG(Value), MAX(Value) FROM T GROUP BY Group

T.Score(3) = NaN;              % introduce missing value

ismissing(T)                   % logical table of missing locations
any(ismissing(T), 'all')      % true if any missing

rmmissing(T)                   % remove rows with any NaN/missing
fillmissing(T, 'linear')       % linear interpolation
fillmissing(T, 'constant', 0) % fill with 0
fillmissing(T, 'movmean', 3)  % moving average

Categorical arrays store discrete, finite sets of values efficiently (like enums). They use less memory than strings and enable grouping operations.

grades = categorical({'A','B','A','C','B','A'});
categories(grades)     % {'A','B','C'}
countcats(grades)      % [3; 2; 1]
grades == 'A'          % logical [1 0 1 0 0 1]

% Ordered categorical:
sizes = categorical({'S','M','L'}, {'S','M','L'}, 'Ordinal', true);
sizes(1) < sizes(2)   % true (S < M)

% Stack (LIFO) using a cell array:
stack = {};
stack{end+1} = 'push item';    % push
item = stack{end};              % peek
stack(end) = [];               % pop

% Queue (FIFO) using a cell array:
queue = {};
queue{end+1} = 'enqueue';    % enqueue
item = queue{1};              % peek front
queue(1) = [];                % dequeue

a = [3 1 4 1 5 9 2 6];
b = [2 7 1 8];

unique(a)             % [1 2 3 4 5 6 9] — sorted unique values
intersect(a, b)       % [1 2] — common elements
union(a, b)           % combined unique elements
setdiff(a, b)         % elements in a but not b
ismember(a, b)        % [0 1 0 1 0 0 1 0] — which a elements are in b

v = 1:2:99;    % sorted odd numbers

% Linear search (O(n)):
find(v == 51)

% Binary search using ismember (optimised internally):
ismember(51, v)      % true

% For insertion point (like Python bisect):
idx = find(v >= 51, 1, 'first')

% histc for binning / bucket search:
edges = [0 25 50 75 100];
histc(v, edges)  % count per bin

d1 = datetime(2024, 6, 15);
d2 = datetime(2024, 12, 31);
d2 - d1                          % duration: 199 days

d1 + days(30)                   % 2024-07-15
month(d1)                        % 6
datestr(d1, 'yyyy-mm-dd')       % '2024-06-15' (char)

% Parse from string:
d = datetime('15/06/2024', 'InputFormat', 'dd/MM/yyyy');

% Array of dates:
dates = datetime(2024,1,1) + calmonths(0:11);

% When output isn't a uniform scalar, use UniformOutput=false:
words = {'hello world', 'foo bar baz', 'matlab'};
parts = cellfun(@(s) strsplit(s), words, 'UniformOutput', false);
parts{1}   % {'hello', 'world'}

% Alternative using a for loop (often clearer and faster):
for k = 1:numel(words)
    parts{k} = strsplit(words{k});
end

% R2019b+: string arrays support many operations directly:
s = ["hello world"; "foo bar"];
split(s)   % splits on whitespace

Fs = 1000;  Fc = 100;    % sample rate, cutoff frequency

% Butterworth lowpass filter (Signal Processing Toolbox):
[b, a] = butter(5, Fc/(Fs/2));   % order 5, normalised cutoff
freqz(b, a, [], Fs)              % frequency response plot

% Apply filter:
x = randn(1, 1000);
y = filter(b, a, x)             % causal (introduces phase delay)
y = filtfilt(b, a, x)          % zero-phase (forward + backward)

% FIR filter via window method:
b_fir = fir1(64, Fc/(Fs/2));   % 64-tap FIR lowpass

Fs = 1000;
t  = 0:1/Fs:1-1/Fs;
x  = sin(2*pi*50*t) + 0.5*randn(1,Fs);

% Welch's method (averaged periodograms):
[pxx, f] = pwelch(x, [], [], [], Fs);
plot(f, 10*log10(pxx));
xlabel('Freq (Hz)'); ylabel('PSD (dB/Hz)');

% Periodogram:
[pxx, f] = periodogram(x, [], [], Fs);

% Continuous wavelet transform (CWT):
[wt, f] = cwt(x, Fs);
surface(t, f, abs(wt));
set(gca, 'YScale', 'log');

% Discrete wavelet transform (DWT):
[cA, cD] = dwt(x, 'db4')      % Daubechies-4 wavelet
xrec = idwt(cA, cD, 'db4')   % reconstruct

% Multi-level decomposition:
[C, L] = wavedec(x, 4, 'db4')  % 4-level

img = imread('cameraman.tif');     % read image (uint8)
imshow(img);                         % display

% Convert to double for arithmetic:
img_d = im2double(img);

% Convert colour to grayscale:
gray = rgb2gray(img);

% Write image:
imwrite(gray, 'gray_output.png');

img = im2double(imread('cameraman.tif'));

% Gaussian blur:
blurred = imgaussfilt(img, 2);

% Edge detection:
edges = edge(img, 'Canny');           % Canny detector
edges = edge(img, 'Sobel');           % Sobel

% Median filter (remove salt-and-pepper noise):
cleaned = medfilt2(img, [3 3]);

% Custom convolution kernel:
sharpen = imfilter(img, fspecial('unsharp'));

bw = imbinarize(imread('coins.png'));

se = strel('disk', 5);       % structuring element
imerode(bw, se)               % erosion (shrink objects)
imdilate(bw, se)              % dilation (expand objects)
imopen(bw, se)                % erode then dilate (remove noise)
imclose(bw, se)               % dilate then erode (fill holes)
bwlabel(bw)                   % label connected components
regionprops(bw, 'Area', 'Centroid')  % measure regions

Fs = 1000;   f0 = 50;   Q = 35;     % quality factor
wo = f0/(Fs/2);
bw = wo/Q;

[b, a] = iirnotch(wo, bw);     % IIR notch filter coefficients
freqz(b, a, 1024, Fs);          % view frequency response

% Apply to noisy signal:
t = 0:1/Fs:1;
x = sin(2*pi*120*t) + sin(2*pi*50*t);  % 120 Hz + 50 Hz hum
y = filter(b, a, x);

Fs = 1000;
t  = 0:1/Fs:2;
x  = chirp(t, 0, 2, Fs/2);    % frequency sweep 0..500 Hz

% stft (R2019a+):
[s, f, t2] = stft(x, Fs, 'Window', hamming(256), ...
    'OverlapLength', 128, 'FFTLength', 512);
imagesc(t2, f, 20*log10(abs(s)));
axis xy; colorbar;
xlabel('Time (s)'); ylabel('Frequency (Hz)');

t = 0:0.01:10;
x = sin(2*pi*0.5*t) + 0.3*randn(1, length(t));

[pks, locs] = findpeaks(x)
[pks, locs] = findpeaks(x, 'MinPeakHeight', 0.5, ...
    'MinPeakDistance', 50)   % at least 0.5 height, 50 samples apart

findpeaks(x, t, 'Annotate', 'extents')  % plot with annotations

img = im2double(imread('cameraman.tif'));

F = fft2(img);           % 2D FFT
Fs = fftshift(F);        % shift DC to centre
imshow(log(1 + abs(Fs)), []);  % log-scaled magnitude
title('2D Fourier Spectrum')

% Inverse:
img_back = real(ifft2(F));   % should match original

% Transfer function: H(s) = (s+2) / (s^2 + 3s + 2)
num = [1 2];
den = [1 3 2];
H = tf(num, den)          % Control System Toolbox

pole(H)                   % [-2; -1] — poles
zero(H)                   % -2 — zeros
bode(H)                   % Bode plot
step(H)                   % step response
impulse(H)                % impulse response
nyquist(H)                % Nyquist diagram

% Plant: double integrator 1/s^2
P = tf(1, [1 0 0]);

% Automatic PID tuning:
C = pidtune(P, 'PID')

% Manual PID controller:
Kp = 10; Ki = 5; Kd = 1;
C = pid(Kp, Ki, Kd)

% Closed-loop step response:
T = feedback(C * P, 1);   % unity feedback
step(T)
stepinfo(T)               % rise time, overshoot, settling time

H = tf([1], [1 2 1]);      % continuous-time system
Ts = 0.01;                  % sample time (seconds)

% Discretise with different methods:
Hd = c2d(H, Ts, 'zoh')     % zero-order hold
Hd = c2d(H, Ts, 'tustin')  % bilinear (Tustin) transform
Hd = c2d(H, Ts, 'matched') % matched pole-zero

% Convert back:
Hc = d2c(Hd, 'tustin')

H = tf([1], [1 2 1]);

isstable(H)                  % true if all poles in left half-plane
p = pole(H);
all(real(p) < 0)             % manual check: negative real parts

% Gain and phase margins:
[Gm, Pm, Wcg, Wcp] = margin(H)
fprintf('Gain margin: %.2f dB, Phase margin: %.2f deg\n', ...
    20*log10(Gm), Pm);

% dx/dt = Ax + Bu,  y = Cx + Du
A = [0 1; -2 -3];
B = [0; 1];
C = [1 0];
D = 0;

sys = ss(A, B, C, D)       % state-space model object

eig(A)                     % eigenvalues = poles of the system
tf(sys)                    % convert to transfer function
step(sys); bode(sys);
lsim(sys, u, t)            % simulate with input u over time t

% Open and configure model:
open_system('myModel');
set_param('myModel', 'StopTime', '10');
set_param('myModel/Gain', 'Gain', '5');

% Run simulation:
simOut = sim('myModel');

% Extract logged signals:
t = simOut.tout;
y = simOut.yout{1}.Values.Data;
plot(t, y);

% Batch parameter sweep:
gains = [1 2 5 10];
for k = 1:numel(gains)
    set_param('myModel/Gain', 'Gain', num2str(gains(k)));
    simOut(k) = sim('myModel');
end

P = tf(1, [1 3 2 0]);    % 1/(s(s+1)(s+2))

rlocus(P);                 % plot root locus

% Find gain for desired closed-loop poles:
[K, poles] = rlocfind(P)  % interactive click on locus

% Programmatically, evaluate at a desired pole location:
s_des = -1 + 2*1i;       % desired dominant pole
K = abs(1/polyval(poly(pole(P)), s_des));

% Continuous-time IIR: freqs (s-domain)
[b, a] = butter(4, 1000*2*pi, 's');  % analog prototype
freqs(b, a)

% Discrete-time: freqz (z-domain)
Fs = 8000;
[b, a] = butter(6, 1000/(Fs/2));    % digital filter
[H, f] = freqz(b, a, 1024, Fs);

subplot(2,1,1); plot(f, 20*log10(abs(H))); ylabel('Mag (dB)');
subplot(2,1,2); plot(f, angle(H)*180/pi);  ylabel('Phase (deg)');

• Growing arrays in loops — pre-allocate with zeros, cell
• Using loops instead of vectorised operations — use element-wise operators
• Calling functions in tight loops — hoist invariant calculations outside
• Unnecessary data copies — avoid modifying large arrays inside functions
• Not using logical indexing — avoids slower find + integer indexing
• Storing data as cell instead of numeric arrays — numerics are faster

% Profile first, optimise second:
profile on; myFunc(); profile viewer

• Use descriptive variable names (signalPower not sp)
• Add a function header comment with description, inputs, and outputs
• Group related code into local functions
• Use constants with UPPER_CASE names
• Avoid magic numbers — assign them to named variables
• Keep functions short and single-purpose

function power = calcSignalPower(signal, sampleRate)
% CALCSIGNALPOWER  Compute RMS power of a signal.
%   power = calcSignalPower(signal, sampleRate)
%   signal     - row vector of samples
%   sampleRate - samples per second (Hz)
    WINDOW_SIZE = round(sampleRate * 0.02);  % 20 ms window
    power = rms(signal);
end

function result = cachedFib(n)
% Fibonacci with memoization using persistent variable
    persistent cache;
    if isempty(cache)
        cache = containers.Map('KeyType','int32','ValueType','double');
    end
    if isKey(cache, n)
        result = cache(n);
        return;
    end
    if n <= 1
        result = n;
    else
        result = cachedFib(n-1) + cachedFib(n-2);
    end
    cache(n) = result;
end

function y = myPlot(x, data, varargin)
    p = inputParser();
    addRequired(p, 'x');
    addRequired(p, 'data');
    addParameter(p, 'Color', 'blue', @ischar);
    addParameter(p, 'LineWidth', 1, @isnumeric);
    addParameter(p, 'Title', '', @ischar);
    parse(p, x, data, varargin{:});
    R = p.Results;

    plot(R.x, R.data, 'Color', R.Color, 'LineWidth', R.LineWidth);
    if ~isempty(R.Title), title(R.Title); end
end

myPlot(x, y, 'Color', 'red', 'LineWidth', 2, 'Title', 'My Plot');

classdef EventBus < handle
    events
        ValueChanged
    end
    properties
        Value = 0
    end
    methods
        function set.Value(obj, val)
            obj.Value = val;
            notify(obj, 'ValueChanged');
        end
    end
end

bus = EventBus;
addlistener(bus, 'ValueChanged', ...
    @(src,~) fprintf('Value changed to: %g\n', src.Value));
bus.Value = 42;   % triggers: "Value changed to: 42"

function y = myFunc(x)
% Works on any size input — uses .* ./ .^ element-wise
    y = sin(x).^2 + cos(x).^2;   % should always be 1
end

myFunc(0.5)            % scalar
myFunc([0 1 2])       % row vector
myFunc(rand(3))        % matrix — all work correctly

function result = divide(a, b)
    if ~isnumeric(a) || ~isnumeric(b)
        error('myPkg:invalidInput', ...
              'Inputs must be numeric; got %s and %s', class(a), class(b));
    end
    if any(b(:) == 0)
        warning('myPkg:divByZero', 'Division by zero detected; returning Inf');
    end
    result = a ./ b;
end

% Structured error ID: 'package:specificError'
% Caller can catch selectively:
try
    divide(1, 0);
catch e
    if strcmp(e.identifier, 'myPkg:invalidInput')
        disp('type error');
    end
end

MATLAB Editor integrates with Git natively (Home → Source Control), or you can use the system command from the Command Window.

% Run any git command from MATLAB:
system('git status')
system('git add myScript.m')
system('git commit -m "Fix indexing bug"')

% Best practices for MATLAB projects in Git:
% - Add .gitignore: *.asv, *.m~ (autosave), *.mex*, slprj/
% - Store figures as PNG/PDF, not .fig files
% - Keep data files separate from code
% - Use .mlx Live Scripts only if needed (binary diff)

Packages are directories starting with +. Functions inside are called with the pkg.func syntax, avoiding name clashes.

% Directory structure:
%  +myPkg/
%    utils.m
%    stats.m
%  main.m

% Inside +myPkg/stats.m:
function m = mean_trim(x, pct)
    lo = prctile(x, pct);
    hi = prctile(x, 100-pct);
    m  = mean(x(x >= lo & x <= hi));
end

% Usage in main.m:
result = myPkg.stats.mean_trim(data, 5);

% Or import:
import myPkg.stats.*
result = mean_trim(data, 5);

% Create a Live Script from script:
matlab.internal.liveeditor.openAndConvert('myScript.m', 'myScript.mlx')

% Export Live Script to PDF or HTML:
% Home → Export → Export to PDF / HTML