%load_ext octavemagic

%octave addpath(genpath("/home/matthew/workspace/resources/N/Numerical Methods/BEARDAHCC/work"))

%%octave

L = [1 0 0 0;2 1 0 0;3 2 1 0;-1 0 2 1];
b = [1 1 1 1]';
L\b

%%octave

U = [-2 1 1 -1;0 -1 0 1;0 0 1 1;0 0 0 2];
y = [1 -1 0 2]';
x = U\y

%%octave

A = [-2 1 1 -1;-4 1 2 -1;-6 1 4 0;2 -1 1 5];
[L,U] = lufac(A,1)

%octave help lufac

%%octave

A = [-2 1 1 -1;-4 1 2 -1;-6 1 4 0;2 -1 1 5];
b = [1 1 1 1]';
[L,U] = lufac(A, 1)
y = L\b
x = U\y

%%octave

newton('newtp5',1,30,1e-3,2,0);

%%octave

type newtp5

%%octave

newton('newtq4',2+2i,30,1e-3,2,0);

%%octave

A7 = [3 0 1;0 -3 -0;1 0 3]
i = 0
x = [1;1;1;]
i = i + 1, x = A7*x

%%octave

i = i + 1, x = A7*x
i = i + 1, x = A7*x
i = i + 1, x = A7*x
i = i + 1, x = A7*x


%%octave

A7 = [3 0 1;0 -3 -0;1 0 3]
i = 0
x = [1;1;1;]
i = i + 1, y = A7*x, [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

i = i + 1, y = A7*x, [z,j]=max(abs(y)), delta=y(j), x=y/delta
i = i + 1, y = A7*x, [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

load matrices
x0 = rand(3,1)
[delta,x,n]=scalepm(M3, x0, 200, 0.001);

%%octave

eig(M3)

%%octave

A7 = [3 0 1;0 -3 0;1 0 3]

%%octave

B = inv(A7)

%%octave

i = 0
x = [0;1;2]
i=i+1, y=B*x, [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

i=i+1, y=B*x, [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

i=i+1, y=B*x, [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

i=i+1, y=B*x, [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

i=i+1, y=B*x, [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

i=i+1, y=B*x, [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

A7 = [3 0 1;0 -3 0;1 0 3]
i = 0
x = [0;1;2]
i=i+1; y=A7\x, [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

i=i+1; y=A7\x, [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

i=i+1; y=A7\x, [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

i=i+1; y=A7\x, [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

i=i+1; y=A7\x, [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

i=i+1; y=A7\x, [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

i=i+1; y=A7\x, [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

i=i+1; y=A7\x, [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

i=i+1; y=A7\x, [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

A7 =[3 0 1;0 -3 0;1 0 3]
[L,U]=lufac(A7,1)
i=0
x=[0;1;2]
i=i+1; y=U\(L\x), [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

i=i+1; y=U\(L\x), [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

i=i+1; y=U\(L\x), [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

i=i+1; y=U\(L\x), [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

i=i+1; y=U\(L\x), [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

i=i+1; y=U\(L\x), [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

i=i+1; y=U\(L\x), [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

i=i+1; y=U\(L\x), [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

A = rand(3)

sigma=1
[V,D]=eig(A)
eig(A)-sigma
[V,D]=eig(A-sigma*eye(3))

%%octave

A7=[3 0 1;0 -3 0;1 0 3]
sigma=-2.5;
M=A7-sigma*eye(3)
i=0
x=rand(3,1)


%%octave

i=i+1, y=M\x, [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

i=i+1, y=M\x, [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

i=i+1, y=M\x, [z,j]=max(abs(y)), delta=y(j), x=y/delta

%%octave

M2=[4 -1 -1 -1;-1 4 -1 -1;-1 -1 4 -1;-1 -1 -1 4]
x0=rand(4,1)
[lambda,x,n]=shiftip(M2,2.9,x0,100,1e-3);

%%octave

lambda

%%octave

circle([4 0],2)

%%octave

hold on
plot(3.9958, -0.2678, 'x')

%%octave

M3 = [4 i -1;2 -5 i;1 0 4+2i]
x0=[1;1;1]
[lambda,x,n]=shiftip(M3,-5,x0,10,1e-3)