p=1;
q=2;
y1=1;
y2=-2;
y=[y1 y2]';
t=(0:0.1:3);

%gaussisk kovariansfunktion
%K='sigma2*exp(-(d/A)^2)';
%eksponentiel kovariansfunktion
K='sigma2*exp(-abs(d)/A)';

%beregn C
A=0.5;
sigma2=3;
C=zeros(2,2);
d=0;
C(1,1)=eval(K);
d=p-q;
C(1,2)=eval(K);
d=q-p;
C(2,1)=eval(K);
d=0;
C(2,2)=eval(K);

%beregn C^-1 y
a=inv(C)*[y1 y2]';

%beregn prediktioner
yhat=zeros(length(t));
r=zeros(1,2);
for i=1:length(t);
 d=t(i)-p;
 r(:,1)=eval(K);
 d=t(i)-q;
 r(:,2)=eval(K); 
 yhat(i)=r*a;
end;

%plotter resultat
plot(t,yhat,'-',[1 2],[1 -2],'+');
