a <- 0
b <- 1
eps <- 0.001

H <-  function(x,y){
  return (- 0.1 * sin (x * (0.5 + y * y)))
}

f <- function(x){
  return (x + 0.1)
}

Integ <- function(x,res,num){
  A <- numeric(length(x))
  A[c(1,n)] <- (x[2]-x[1])/2
  A[2:(n - 1)] <- (x[2]-x[1])
  result <- numeric(1)
  for(i in 1:length(x))
    result <- result + A[i]*H(res[num],x[i])*res[i]
  return (result)
}

n <- 5
repeat {
  h <-  (b - a)/(n - 1)

  x <- numeric(n)
  for(i in 1:n)
    x[i] <- a + h * (i - 1)

  f_mat <-  matrix(f(x))

  H_mat <- matrix(nrow = n, ncol = n)
  for(i in 1:n)
    for(j in 1:n)
      H_mat[i,j] <-  H(x[i],x[j])

  A <- numeric(n)
  A[1] <- h/2
  A[n] <- h/2
  A[2:(n - 1)] <- h

  D <- matrix(nrow = n, ncol = n)
  for(i in 1:n)
    for(j in 1:n)
      D[i,j] <- ifelse(i == j, 1, 0) - A[j] * H(x[i],x[j])

  res <- solve(D, f_mat)
  finish <- T
  print(n)
  for(i in 1:3){
    num <- 1 +(i - 1)*(n - 1)/2
    finish <- finish & (abs(res[num] - Integ(x, res, num)-f(res[num])) < eps)
    print(res[num])
    print(abs(res[num] - Integ(x, res, num)-f(res[num])))
  }
  if(finish)
    break
  n <- (n - 1) * 2 + 1
}