Sistemi lineari con metodo delle tangenti matlab

Ciao a tutti, sono nuovo e visto che sul web non c'è molto ho deciso di rivolgermi a voi
Da qualche giorno sono bloccato sul seguente problema:
ho un sistema lineare di 4 equazioni in 4 incognite, e fin qui tutto semplice se non fosse che le incognite sono angoli argomenti di funzioni trigonometriche, per la risoluzione quindi ho usato il metodo iterativo delle tangenti o di newton-raphson.
Sono riuscito a scrivere un programma che mi da le soluzioni (che chiamerò B e C) delle prime due equazioni e queste sono corrette. Il problema sta nel fatto che le soluzioni vengono fuori in base ad un valore che io assegno (un altro angolo chiamato A), dunque per ogni valore di A mi vengono fuori un valore di B ed uno di C; tuttavia non riesco a dire al programma di prendere e leggere i valori di A (che varia tra 0 e 360) e per ognuno darmi il valore di B e C in modo da ottenere in output un grafico di B e C in funzione di A (che è lo scopo di tutto).
Ho strutturato il programma in quattro moduli:
1.fun dove dichiaro le variabili e le equazioni
2.jacobian dove scrivo lo jacobiano associato al sistema
3.il metodo di newton raphson
4.il main (che secondo me è quello in cui c'è il problema)
Potete aiutarmi?
Grazie mille a tutti!

modulo fun

function f=fun(x)
% 
% INPUT
% x: Vettore delle incognite
% OUTPUT
% f: Vettore dei residui

r(2)=1;
r(3)=3;
r(4)=4;
r(6)=4;
r(7)=3;
for i=1:1:360	%qui ho sicuramente sbagliato qualcosa perchè in output stampa i 360 valori ti theta
    theta(1)= i;
end

f(1)= r(2)*cosd( theta(1)) + r(3)*cosd( x(1)) + r(4)*cosd( x(2)) - r(6);
f(2)= r(2)*sind( theta(1)) + r(3)*sind( x(1)) + r(4)*sind( x(2)) + r(7);
f= f';

end

modulo jacobiano

function J = jacobian(x)
%
% INPUT
% x: Vettore delle incognite
% 
% OUTPUT
% J: Matrice Jacobiana
r(3)=3;
r(4)=4;
r(2)=1;

J(1,1) = -r(3)*sind( x(1));
J(1,2) = -r(4)*sind( x(2));
J(2,1) =  r(3)*cosd( x(1));
J(2,2) =  r(4)*cosd( x(2));

end

modulo newton raphson

function [x,iterations,f,resid,exitflag]=newtonraphson(fun,jacobian,x0,tol,itmax)
%
% INPUT
% fun : Function in cui vengono valutati gli elementi del vettore dei
% residui(equazioni di chiusura)
% jacobian : Function in cui vengono valutati gli elementi della matrice jacobiana
% x0: Vettore soluzione di primo tentativo(i valori misurati che io do in input)
% tol: Tolleranza (precisione con cui voglio la soluzione)
% itmax: Numero massimo di iterazioni
%
% OUTPUT
% x: Soluzione
% iterations: Numero di iterazioni eseguite
% f: Vettore dei reidui (equazioni di chiusura)
% resid: Norma del vettore dei residui
% exitflag:  1 L'iterazione converge sul valore x, 0
%            0 Si e' superato il massimo numero di iterazioni

% Inizializza le variabili
resid=1000;
iterations=0;
while(resid>tol)
    % Calcola matrice Jacobiana
    J=feval(jacobian,x0);
    % Calcola residui
    f=feval(fun,x0);
    % Calcola norma del vettore dei residui
    resid=norm(f);
    
    % Trova il nuovo valore di x con il metodo di Newton-Raphson
    x = x0-J\f; % inv(J)  J \ 
    
    
    % Controllo sulla tolleranza
    if(resid<tol)
        exitflag=1;
        return
    end
    % Controllo sul numero delle iterazioni
    if(iterations>itmax)
        exitflag=0;
        return
    end
    % Assegan il nuovo valore come condizioni iniziali per la successiva
    % iterazione
    x0=x;
    iterations=iterations+1;
end

ed infine il main

clc; clear; close all;
x0=[257.85 ; 348.84]; 


str = sprintf('Soluzione di primo tentativo: %.1f  %.1f',x0);
disp(str);

% Criterio convergenza
tol=1e-16;
str = sprintf('Tolleranza: %.2e',tol);
disp(str);

% Massimo numero di iterazioni
itmax=100000;
str = sprintf('Massimo numero di iterazioni: %i',itmax);
disp(str);


str= sprintf ( '\n\n----Metodo Newton-Raphson----\n\n');
disp(str);
% Metodo di Newton-Raphson per trovare la soluzione del sistema di
% equazioni non lineari
[x,iterazioni,funzioni,residui,exitflag]=newtonraphson(@fun,@jacobian,x0,tol,itmax);


str = sprintf('Soluzione = %.2e %.2e',x);
disp(str);
str = sprintf('Iterazioni = %.i',iterazioni);
disp(str);
str = sprintf('Residui della funzione = %.3e %.3e',funzioni);
disp(str);
str = sprintf('Massima norma dei residui = %.3e',residui);
disp(str);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%                    PLOT                              %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r(2)=1;
r(3)=3;
r(4)=4;
r(6)=4;
r(7)=3;
%theta(2)=pi/4;
%theta(4)= pi;
for i=1:1:360	%qui come prima non sono sicuro sia giusto
    theta(1)=i
end
f(1)= r(2)*cosd( theta(1)) + r(3)*cosd( x(1)) + r(4)*cosd( x(2)) - r(6);
f(2)= r(2)*sind( theta(1)) + r(3)*sind( x(1)) + r(4)*sind( x(2)) + r(7);
z_point= r(2)*cosd( theta(1)) + r(3)*cosd( x(1)) + r(4)*cosd(x(2)) + r(6);




figure('Position',[100 100 900 500])
colormap jet % jet copper cool
title('\fontsize{15}Grafico delle posizioni')
grid on
xlabel('Theta1','fontweight','b');
ylabel('theta3 e theta2','fontweight','b');
hold on
plot( theta(1), x(1),'ro','Linewidth',4)
plot(theta(1), x(2), 'ro','Linewidth',4)