Fortran to Matlab, B-Spline,

program bspline_fit

    use toolbox

    implicit none
    integer :: ord,typ,ni,no
    real :: t,dt
    real, allocatable :: nvo(:),nvi(:),crd0(:,:),crd1(:,:)
    integer :: i,j

    ! input data (the number of control points will be equal to the number of data points)
    write (*,*) "Enter the number of data points"
    read (*,*) ni

    write (*,*) "Enter the order of the curve"
    read (*,*) ord
    no = ni+ord

    write (*,*) "Enter the type of the knot vector"
    write (*,*) "(1) periodic uniform (2) open uniform (3) nonuniform"
    read (*,*) typ

    allocate (crd0(ni,2),crd1(ni,2),nvi(ni),nvo(no))

    call data_read(ni,crd0) ! input data
    call data_vec(ni,crd0,nvi) ! parametric coordinates of the data points
    call knot_vec(typ,ord,ni,no,nvo) ! knot vector
    call control_vec(ord,ni,no,nvi,nvo,crd0,crd1) ! parametric coordinates of the control points
    call curve(ord,ni,no,nvi,nvo,crd1) ! resulting curve

end program bspline_fit

--------------------------------------

module toolbox

    contains

    subroutine data_read(n,crd)

    implicit none
    integer,intent(in) :: n
    real,intent(out) :: crd(n,2)
    integer :: i,j,k

    open (1,file = 'in.dat')
    do i=1,n ; read (1,*) k,(crd(i,j),j=1,2) ; end do

    end subroutine data_read

    subroutine data_vec(n,crd,nvi)

    implicit none
    integer,intent(in) :: n
    real,intent(in) :: crd(n,2)
    real,intent(out) :: nvi(n)
    real :: dif(2),clen(n-1)
    integer :: i

    do i=1,n-1
        dif = crd(i+1,:)-crd(i,:)
        clen(i) = sqrt(dif(1)**2+dif(2)**2)
    end do
    nvi(1) = 0
    do i=2,n
        nvi(i) = nvi(i-1)+clen(i-1)
    end do
    nvi = nvi/sum(clen)
    write (*,'(20(f5.2,2x))') (nvi(i),i=1,n) ; write (*,*)

    end subroutine data_vec

    subroutine knot_vec(typ,ord,ni,no,nvo)

    integer,intent(in) :: typ,ord,ni,no
    real,intent(out) :: nvo(no)
    integer :: i

    select case(typ)
    case(1)

    case(2)
        do i=1,no
            if (i <= ord) then
                nvo(i) = 0
            else if (i >= ni+1) then
                nvo(i) = ni+1-ord
            else
                nvo(i) = i-ord
            end if
        end do
    case(3)

    end select
    nvo = nvo/nvo(no)
    write (*,'(20(f5.2,2x))') (nvo(i),i=1,no)

    end subroutine knot_vec

    subroutine control_vec(ord,ni,no,nvi,nvo,crd0,crd1)

    integer,intent(in) :: ord,ni,no
    real,intent(in) :: nvi(ni),nvo(no),crd0(ni,2)
    real,intent(out) :: crd1(ni,2)
    real :: cmat(ni,ni)
    integer :: i

    open (3,file='tmp.dat')
    do i=1,ni
        call basis_fnc(ord,ni,no,nvo,nvi(i),cmat(i,:))
    end do
    do i=1,ni ; write (*,'(20(f6.3,2x))') (cmat(i,j),j=1,ni) ; end do ; write (3,*)
    call mat_inv(ni,cmat) ; crd1 = matmul(cmat,crd0)
    do i=1,ni ; write (3,'(i3,1x,3(f6.3,2x))') i,crd1(i,1),crd1(i,2) ; end do ; write (3,*)

    end subroutine control_vec

    subroutine basis_fnc(ord,ni,no,x,t,cv)

    implicit none
    integer,intent(in) :: ord,ni,no
    real,intent(in) :: x(no),t
    real,intent(out) :: cv(ni)
    real :: bf(ord,no-1),r
    integer :: i,j,k

    bf = 0
    k = no-1
    do i=1,ord
        if (i == 1) then
            do j=1,k
                if (t >= x(j) .and. t <= x(j+1)) then
                    bf(i,j) = 1
                else
                    bf(i,j) = 0
                end if
            end do
        else
            k = k-1
            do j=1,k
                r = x(i+j-1)-x(j) ; if (r /= 0) bf(i,j) = (t-x(j))*bf(i-1,j)/r
                r = x(i+j)-x(j+1) ; if (r /= 0) bf(i,j) = bf(i,j)+(x(i+j)-t)*bf(i-1,j+1)/r
            end do
        end if
        write (3,'(20(f6.3,2x))') (bf(i,j),j=1,k)
    end do
    cv = bf(ord,1:ni)

    end subroutine basis_fnc

    subroutine curve(ord,ni,no,nvi,nvo,crd1)

    implicit none
    integer,intent(in) :: ord,ni,no
    real,intent(in) :: nvi(ni),nvo(no),crd1(ni,2)
    real :: cv(ni),crd(2),t
    integer :: i,j,k
    real,parameter :: dt = 1./20

    open (2,file = 'out.dat')
    t = 0 ; i = 1
    do
        call basis_fnc(ord,ni,no,nvo,t,cv)
        do j=1,2 ; crd(j) = dot_product(cv,crd1(:,j)) ; end do
        write (2,'(i3,2x,2(f5.2,2x))') i,(crd(j),j=1,2)
        i = i+1 ; t = t+dt ; if (t > nvi(ni)) exit
    end do

    end subroutine curve

    ! matrix inverse
    subroutine mat_inv(n,a)

    implicit none
    integer,intent(in) :: n
    real,intent(inout) :: a(n,n)
    integer :: i
    real :: v1(n),v2(n)

    call sgetrf(n,n,a,n,v1,i) ! lapack lu decomposition
    call sgetri(n,a,n,v1,v2,n,i) ! lapack matrix inverse

    end subroutine mat_inv

end module toolbox

---------------------------------------

Is there anyone can help me to transform this code to Matlab?