function [A,b] = mkd1d_NeumannVpravo(p,f,tau,l,n)

h = l/(n+1);
B = ones(n+1,1)*[-1,2,-1];
A = (p/h^2) * spdiags(B,-1:1,n+1,n+1);
A(n+1,n:(n+1)) = (p/h^2)*[-h,h];
x = h*[1:n]';
b = [f(x);tau];